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ABSTRACT 

The halos of galaxies preserve unique records of their formation histories. We carry out the first 
combined observational and theoretical study of phase-space halo substructure in an early-type galaxy: 
M87, the central galaxy in the Virgo cluster. We analyze an unprecedented wide-field, high-precision 
photometric and spectroscopic data set for 488 globular clusters (GCs), which includes new, large- 
radius Subaru/Suprime-Cam and Keck/DEIMOS observations. We find signatures of two substruc- 
tures in position-velocity phase-space. One is a small, cold stream associated with a known stellar 
filament in the outer halo; the other is a large shell-like pattern in the inner halo that implies a mas- 
sive, hitherto unrecognized accretion event. We perform extensive statistical tests and independent 
metallicity analyses to verify the presence and characterize the properties of these features, and to 
provide more general methodologies for future extragalactic studies of phase-space substructure. The 
cold outer stream is consistent with a dwarf galaxy accretion event, while for the inner shell there 
is tension between a low progenitor mass implied by the cold velocity dispersion, and a high mass 
from the large number of GCs, which might be resolved by a ~ 0.5 i* E/SO progenitor. We also 
carry out proof-of-principle numerical simulations of the accretion of smaller galaxies in an M87-like 
gravitational potential. These produce analogous features to the observed substructures, which should 
have observable lifetimes of ~ 1 Gyr. The shell and stream GCs together support a scenario where 
the extended stellar envelope of M87 has been built up by a steady rain of material that continues 
until the present day. This phase-space method demonstrates unique potential for detailed tests of 
galaxy formation beyond the Local Group. 

Subject headings: Galaxies: ellipticals and lenticulars, cD — Galaxies: halos 



1. INTRODUCTION 

The ultimate record of the formational histories of 
galaxies is found in the six-dimensional phase space of 
stellar positions and velocities, supplemented by addi- 
tional information about ages and elemental abundances 
(primarily iron-based metallicities) . The halos of galaxies 
are natural hunting grounds for assembly clues in phase- 
space because of the preservative effects of the long dy- 
namical times, and since any material entering a galaxy 
naturally has to traverse its halo. 

Indeed, it is a fundamental prediction of the current 
cosmological paradigm that the outer regions of galax- 
ies should be constantly bombarded with smaller, in- 
falling systems that persist as significant substructures 
while gradually becoming disrupted. This idea has been 
borne out qualitatively in recent years by observations 
of halo streams and shells in the nearb y Universe (e.g., 
Helmi et all 119 99': 'Ib ata et al.l '2001: Bclok urov et alJ 
20061: iTal eral., ,2009^ iMartinez-Delgado et all 120101 : 
Cooper et al.!i20m 

Massive elliptical galaxies are touchstones for this is- 
sue, since their spheroidal nature and their residence 
in the densest areas of the Universe suggest partic- 
ularly vigorous accretion histories. Comparisons of 
these galaxies' stellar luminosity profiles at high and 



low redshifts now provide indirect evidence that their 
outer envelopes grow tremendously in size up to the 
present dav fe.g.jBu itrago ct al." 2008; van der Wei et al] 
20081 iDamianov et al., .2009.; van Dokkum et al.l 120101 
Cassata et al.ll2011[ ). which may be naturally explained 
by accretion and mergers in a c osmologic al context 
jKhochfar & Silk 2006; Naa b et all [2009: Hop kins et all 
[20101: lOser et al.li2010l . l201lD . 

The most extreme examples of ongoing assembly 
are expected to be the brightest cluster galaxies 
(BCGs), whose typical locations at the centers of clus- 
ters should involve very ac tive merger histories (e.g., 
iRuszkowsk i fc Springelll2009f )r'1 Current cosmologically- 
bascd models find that around half of the stellar mass 
growth of a present-day BCG occurred over the past 
5 Gyr via gas-poor merg ers (after an early forma tion 
epoch for the core regions; IDe Lucia fc Blaizot|[2007l) . 

This theoretical picture has received very mixed sup- 
port from comparisons of BCG stellar masses and 
surface brightness distributions at diff erent redshifts 
(iWhilev et al.|[20q8l: [Bernardi 2009; CoUins et al.ll2009al: 
iValentinuzzi et al.l 120101 : lAscaso et al.l 120111 : iStott et al.l 

^ Although the central galaxy is often not technically the bright- 
est object in a cluster, we will follow the convention of calling this 
the BCG. 
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1201 It) . Surveys for niajor mergers in action have 
found th at these contribute si g nificantly to rec ent BCG 
growth (iMcIntosh et all 120081: iLiu et all l2009t see also 
IBroughet all 120 111) , but the contributions from minor 
mergers are unknown and the overall comparison to the- 
ory unclear. Part of the difficulty is that even if accre- 
tion and mergers are important growth processes, the 
frequency of observing an ongoing event can be fairly 
low. Thus it would be invaluable to go beyond purely 
photometric observations and delve into the phase-space 
of position, velocity, and mct allicit y for longer lived sig- 
natures (e.g.. [Johnston et ani200l ' lZemp et al.H200l . 

Within the Local Group, detailed phase-space stud- 
ies are being undertaken using ind i vidual stars (e.g., 
Koch et a l.|[2008t [Gilbert et al.ll200l iStarkenburg et alj 



2009; Xu e et al.ll20lH ). but are currently impossible in 
more distant galaxies. Instead, bright stellar prox- 
ies may be used, such as planetary nebulae (PNe; 
Durrell et ahl [20031 iMeFrett et all [2001 iShih fc Mended 



20101: iCortesi et al.ri201ir and globular clusters (GCs). 
GCs are thought to form along with field stars, and 
as collections of ~ 10^ stars each, they are visible at 
much greater distances, allowing their individual posi- 
tions, line-of-sight velocities, and chemica l properties to 
be m easured as far away as ^ 50 Mpc (|Richtler et al.l 

[mil). 

The "chemodynamics" of GCs provided the water- 
shed evidence for the accretion origin of the Milky 
Way's outer halo (jSearle &: ZinnI [1978!), where it is 
now thought that disrupting satellite galaxies h ave de- 
posited many accompanying GCs (e.g., [Bcllazzin i et al.l 



i ccompanymg WL-s (e.g., |B 

200llG ao et al."2007( ' Geisler et all^OOTtTLee et alll2007 



Forbes & Bridges 2010). A similar scenario is now ev- 



ident in the nearby spiral M31, where some GCs are 
associated with hal o substructures (Collins et al,, 2009tx 
iMackev et aT1[20Tol ). In BCGs, the presence of GCs in 
enormous numbers (up to 5 x lO'* per system) has long 
motivated suggestion s that accretion is an important fac - 
tor for these galaxies (|Forte et al.lll982l : lC6te et al.lll998f ). 
although it is not clear how many of these GCs may have 
formed during in situ processes or mergers rather than 
being accreted. 

Various GC and PN-based studies of the halo kinemat- 
ics of giant ellipticals and BCGs have turned up evidence 
for substructur es that may ir nply a ct ive accretion events 
or me rgers (iCote et al.l 120031 : iRomanowskv et al.l 



20091: iSchuberth et al 
Woodlev fc Harris! 1201 ID . 



l2010t [McNeil et al.l 12'oTor 

In other cases, accretion 



events are suggested by photometric substructures or by 
broad kinematical or chemi c al halo transitio ns ( Tal ct al. 
200^; 'Proctor et al.' '200^; 'Coccato et al.' '2009,j2010t 
Janowiccki ct al. 2010 ' Forbes et al. 2011; Arno ld et al.l 
20ilt iMcmhc me et al.l l201lT) ! However, in no case has 



a halo feature been observed and modeled in enough 
detail to determine its origin and the implications for 
galaxy formation. 

To go beyond the previous work on nearby BCGs by 
carrying out detailed phase-space studies of their ha- 
los, there are two critical observational requirements be- 
side sheer instrumental throughput: wide-field coverage 
and spectroscopic resolution. The former permits blind 
searches for halo substructures and demands a field of 
view of at least ~ 0.5 degree in order to probe galac- 



tocentric radii out to ~ 100 kpc or more in galaxies at 
~ 15 Mpc distances. The latter allows for a wide mass 
range of accreted galaxies to be probed; e.g., if a BCG 
has a characteristic velocity dispersion of cr 300 km s~^ 
and the precision of velocity measurements for halo trac- 
ers is Av ^ 50-100 km s~^, as in many past surveys, 
then the observations are sensitive to mass ratios down 
to ~ (Aw/cr)^ ~ 1:40-1:10 in dynamical mass. More 
minor mergers than these are certainly expected to be 
more frequent, and may comprise a dominant mode for 
halo assembly, so a velocity resolution of tens of km s~^ 
(Aw/cr ^ 0.1) is preferable. 

We present here the first wide-field, large-sample spec- 
troscopic survey of a BCG with velocity resolution of 
Av/a 0.04. Our subject is M87, the central el- 
liptical in Virgo, the nearest galaxy cluster at a dis- 
tance of ^ 16.5 Mpc. The GC system of M87 ha s seen 
decades of spectroscopic studv fe.g.. iHuchra fc Brodid 
[19871: [Cohe"n fc Rvzhovl[T997l: IHanes et al.ll200lD . but ex- 
tending to radii of only ~ 40 kpc and with Av/a ^ 0.25. 

Until now, M87 could be characterized as a dynami- 
cally quiet galaxy, with only mild signs of interactions 
such as very faint s tellar filaments in its fa r outer halo 
(jMihos et al.l [20051 : [jaM)wiecki et al.l IMol ). O ur new 
GC study reveals the kinematics of one of these faint 
streams, and unveils an unsuspected, enormous shell-like 
substructure in phase-space. 

Our paper proceeds as follows. Section [2] describes the 
observations and data reduction, and presents a basic 
overview of the halo substructures. Sections [3| and [4| an- 
alyze the characteristics and statistical significance of the 
large inner shell and outer stream, respectively. Theoret- 
ical analyses of the substructure origins and dynamics are 
explored in Section [5l Section [6| summarizes the findings 
and outstanding questions. 

2. OBSERVATIONS AND BASIC RESULTS 

We have revisited M87 in a new-era high-precision 
wide-field survey of GCs. The main component of 
our data set i s pre sented and discussed in detail in 
iStrader et all (|2Cilll . hereafter S+11). In brief, it 
is based on photometry from CFHT/Megacam and 
Subaru/Suprime-Cam optical images, and spectroscopy 
from MMT/Hectospec, Keck/DEIMOS, and Keck/LRIS. 
New line-of-sight velocities were obtained for 451 GCs 
over a series of campaigns from 2007 to 2010, extending 
over a range in galactocentric radius of ~ 1-35 arcmin 
(~ 5-170 kpc). 

The first subset of the new data extended along a nar- 
row track eastwards of M87's center, in search of any 
kinematical transition between the BCG and the sur- 
rounding cluster. A peculiar fluctuation in the velocity 
dispersion proflle was indeed found at a radius of ^ 9 ar- 
cmin (~ 45 kpc), which was the first indication of an 
inner-halo shell-like feature as we will discuss later. Sub- 
sequent data (the large majority of the total sample) were 
obtained at generally random position angles around the 
galaxy. 

In a parallel campaign presented here for the first 
time, we observed in tensively the region around the outer 
stellar "stream A" (iMihos et al.l [20051: iJanowiecki et al.l 
[2OT0t IRudick et al]l2010l: IKrick et al.ll201lD . at radn of 
~ 35-50 arcmin (~ 170-240 kpc). GC selection for spec- 
troscopic follow-up was obtained from a variety of images 
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as they became availab le: archival Subaru/Suprime-Cam 
(jMivazaki et al.|[200l BVI imaging, th e Megacam gri 
imaging mentioned above (jHarrisI I2009D , and our own 
new Suprime-Cam imaging in gi. Given the incomplete 
spatial coverage of the first two imaging data-sets, the 
third data-set provides our default photometric source 
for the stream area. These images were taken on 2009 
Apr 21, with intermittent transparency, 1" seeing, and 
exposure times of 300 sec in each band. 

Our spectroscopic follow-up of the stream GCs 
used Keck/DEIMOS and observing and analysis tech- 
niques e stablished in our previou s work on GC kine- 
matics (iRomanowskv et al.l 120091: lArnold et all 120111: 
[Foster et al.l 120111 : S-t-ll) Four masks were observed 
during the nights 2009 Mar 23, 2010 Mar 11-13, and 
2010 Jun 12, with generally good conditions and expo- 
sure times varying from 0.5 to 2 hours. The spectra 
cover an approximate spectral range of 6500-9000 A at 
a moderate resolution (1.5 A FWHM) that allows the 
derivation of precise line-of-sight velocities. 

The data were reduced in a standard manner, includ- 
ing bias subtraction, flat fielding, wavelength calibra- 
tion, and sky subtraction, followed by optimal extraction. 
Heliocentric line-of-sight velocities were derived through 
cross-correlation with a range of stellar templates, us- 
ing only the region around the Calcium II triplet (CaT; 
Ha, if available, was used for confirmation of borderline 
cases) . Some background galaxies were identified via ex- 
amination of the extracted spectra and excluded from 
further analysis. 

The velocity uncertainties were estimated by com- 
bining in quadrature the statistical uncertainties from 
the cross-correlation with the systematic uncertainties 
from using different stellar templates. The uncertain- 
ties for the confirmed GCs ranged from 5 to 18 km s~^, 
with their reliability supported by repeat observations 
of four GCs and five stars. The characteristics of the 
spectroscopically-confirmed GCs (20 of them), stars, and 
galaxies in the stream region are provide d in Table [U We 
also t argeted five PN candidates from iFeldmeier et al.l 
(|2003| ) but found no Ha emission in their spectra (ob- 
jects 7-14, 7-18, 7-51, 7-58, 7-65). 

TABLE 1 

Spectroscopic ob.iects in M87 outer stream region 
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ID 


R.A. 


DecL 




(9- 


- i)o 


V 




[J2000] 


[J2000] 








[ km s"-'] 


confirmed: 
















R8922 


187.36271 


12 


84842 


20.07 


1 


01 


1460 ± 6 


R10106 


187.35946 


12 


87679 


19.52 





99 


1500 ± 5 


R10806 


187.36816 


12 


89351 


21.36 


1 


08 


839 ±8 


R10814 


187.35690 


12 


89390 


21.38 





85 


1386 ± 10 


R11987 


187.36719 


12 


92166 


22.24 





78 


1100 ±10 


R13629 


187.32990 


12 


96460 


21.65 





80 


1112 ± 7 


R14045 


187.35667 


12 


97587 


21.72 





80 


852 ±8 


R14100 


187.29791 


12 


97751 


21.77 





97 


1154 ±11 


R15088 


187.32974 


13 


00221 


21.09 





76 


1395 ± 10 


R15675 


187.27000 


13 


01773 


21.98 





84 


1141 ± 18 


R16138 


187.41290 


13 


02971 


22.07 


1 


13 


1063 ± 12 


R16501 


187.36578 


13 


03980 


22.33 





71 


1590 ±12 


R16571 


187.36916 


13 


04153 


20.44 





79 


1577 ± 9 


R17084 


187.33063 


13 


05539 


21.18 





83 


1490 ± 8 


R17186 


187.23102 


13 


05824 


20.41 





84 


1291 ± 6 


R17522 


187.34823 


13 


06668 


20.89 





78 


1467 ± 5 


R18039 


187.34714 


13 


08054 


22.08 





76 


673 ±9 


R18045 


187.31874 


13 


08087 


21.47 





81 


1462 ±11 


R19909 


187.33585 


13 


12520 


21.79 





83 


1437 ± 12 



ID 


R.A. 


Dccl. 


io 


(9 - i)o 


V 




r 190001 


[190001 






^ kill s ' ] 


R21956^ 


187.23676 


13.16944 


19.90 


0.87 


1876 ± 7 


marginal: 












R12580 


187.39820 


12.93529 


22.01 


0.74 


1085 ± 16 


R13657 


187.39522 


12.96553 


22.14 


0.78 


1222 ± 18 


R15867 


187.41155 


13.02209 


22.31 


0.75 


1244 ± 14 


R16402 


187.32270 


13.03704 


22.45 


1.21 


898 ± 17 


stars'^ : 












R7403 


187.37441 


12.81474 


17.55 


0.59 


-3 ±48 


R7754 


187.37956 


12.82188 


19.87 


1.91 


14± 74 


R8193 


187.38146 


12.83134 


19.29 


0.69 


-148 ± 21 


R8376 


187.37963 


12.83570 


20.52 


0.98 


-134 ±23 


R8798 


187.39239 


12.84514 


18.86 


0.72 


65 ±21 


R9363 


187.34117 


12.85885 


19.34 


0.88 


40 ±7 


R10146 


187.34439 


12.87770 


20.45 


0.65 


58 ±20 


R10837 


187.36675 


12.89449 


17.67 


2.65 


42± 14 


R11062 


187.40232 


12.89997 


15.38 


1.74 


-24 ±5 


R11805 


187.39709 


12.91762 


22.47 


1.33 


189 ± 20 


R11872 


187.36193 


12.91931 


19.96 


1.11 


45± 15 


R12125 


187.40527 


12.92458 


20.42 


0.74 


145 ± 6 


R12313 


187.29011 


12.92963 


15.75 


0.49 


4 ±20 


R12657 


187.34624 


12.93754 


17.48 


1.09 


10 ±5 


R12691 


187.40031 


12.93805 


16.32 


2.03 


-20 ±5 


R12891 


187.34345 


12.94330 


20.44 


0.84 


177 ±51 


R13164 


187.30174 


12.95084 


20.11 


0.69 


-41 ± 25 


R13188 


187.23932 


12.95184 


17.91 


1.33 


35± 17 


R13506 


187.46370 


12.95930 


15.37 


0.35 


-6 ±8 


R13992 


187.33976 


12.97469 


19.83 


1.80 


4 ±36 


R14054 


187.45815 


12.97566 


19.07 


0.74 


98 ±5 


R14166 


187.33685 


12.97877 


20.88 


1.03 


198 ± 7 


R14317 


187.34125 


12.98226 


15.49 


0.92 


30 ±5 


R14391 


187.35850 


12.98435 


19.70 


0.69 


-87 ±6 


R14427 


187.44737 


12.98488 


20.96 


1.03 


213 ± 14 


R14800 


187.26911 


12.99484 


15.70 


1.67 


0± 18 


R15119 


187.40424 


13.00259 


17.83 


2.02 


9±5 


R15179 


187.44657 


13.00379 


17.83 


0.87 


70 ±5 


R15494 


187.45031 


13.01196 


19.09 


0.82 


-3 ±5 


R15721 


187.30143 


13.01873 


16.86 


0.59 


58± 18 


R15838 


187.42486 


13.02142 


21.09 


0.72 


-31 ± 10 


R15983 


187.36178 


13.02593 


21.37 


1.89 


55 ±32 


R16488 


187.39562 


13.03932 


19.52 


1.12 


108 ± 5 


R17262 


187.41103 


13.05953 


22.24 


0.64 


110 ±24 


R17841 


187.35704 


13.07556 


22.00 


0.72 


185 ± 12 


R18139 


187.27810 


13.08366 


22.26 


1.27 


-84 ± 26 


R18154 


187.34945 


13.08374 


18.48 


1.54 


-57 ±5 


R18730 


187.27131 


13.09786 


21.26 


0.80 


-214 ± 15 


R19052 


187.27947 


13.10543 


15.40 


0.97 


-35 ± 15 


R19452 


187.37246 


13.11475 


16.77 


2.38 


-1±5 


R20070 


187.26784 


13.12938 


17.82 


0.67 


127 ±14 


R20445 


187.30037 


13.13758 


20.15 


0.76 


103 ± 6 


R20496 


187.38952 


13.13803 


21.20 


0.78 


32± 16 


R20765 


187.24172 


13.14418 


20.54 


0.72 


222 ± 18 


R21009 


187.23749 


13.14917 


15.18 


1.85 


5± 17 


R21082 


187.29307 


13.15058 


17.69 


1.30 


-6 ±8 


R22195 


187.33737 


13.17345 


20.94 


0.79 


73 ±9 


R22421 


187.33604 


13.17731 


21.49 


0.57 


143 ± 16 


R22843 


187.32868 


13.18426 


20.83 


0.71 


22 ± 13 


galaxies: 












R12858 


187.23866 


12.94289 


22.24 


0.71 




R13471 


187.25251 


12.95961 


21.89 


1.05 




R14229 


187.22466 


12.98077 


22.03 


0.94 




R15594 


187.25608 


13.01557 


22.41 


0.98 




R15685 


187.23852 


13.01819 


22.40 


1.37 




R17134 


187.23876 


13.05692 


22.59 


0.77 




R21550 


187.21381 


13.15999 


22.97 


0.69 





2 Probably bound to NGC 4461. 
^ A simple velocity boundary of 250 km 



has been used to 

separate stars and GCs, but some of the higher- velocity stars might 
in principle be low- velocity GCs. 
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TABLE 1 — Continued 



ID 


R.A. 


Dccl. 


«o 


(9 - «)o 


V 




[J2000] 


[J2000] 






[km s-i] 



Note. — Photometry is from Subaru-Suprime/Cam imaging 
using 4-pixel aperture magnitudes, bootstrapped to the Sloan Dig- 
ital Sky Survey and then corrected for Galactic extinction follow- 
ing jiJeelc^2iMii {sold). The uncertainties in the stellar veloc- 
ities are indicative only, and not as carefully characterized as for 
the GCs. The most likely objects belonging to the cold GC stream 
are R16501, R16571, R17084, R17522, R18045. 

Our new spectroscopic data thus provides a total of 468 
velocities of M87 GCs (we exclude one large-radius object 
probably associated with the Virgo galaxy NGC 4461, 
but retain another that may be bound to the small el- 
liptical NGC 4478). Although a great deal of additional 
spectroscopic data of M87 GCs are available in the lit- 
erature, we do not use most of the previous data, since 
the typically large measurement errors impede the study 
of cold substructures. Instead, we incorporate only 20 
high resolution measurements from the literature (with 
median uncertainties of 4 km s~^: iHasegan et al.l 120051 : 
lHa^egan|[2007t lEvstigneeva et al.|[2007[) . Note that our 
treatment of previous data differs slightly from the con- 
vention of S-t-11, who incorporated all data published 
since 2003. 

The combined literature and new data provide a total 
sample of 488 GC line-of-sight velocities around M87, ex- 
tending from its central regions out to 200 kpc in its halo, 
where the median velocity uncertainties are 14 km s~^. 
This contrasts with most previous work on GC kinemat- 
ics beyond the Local Group, with typical velocity errors 
of ^ 30-100 km s^^. Our GC sample also has an ex- 
tremely low rate of contamination from foreground stars 
and background galaxiesQ This is the largest GC data- 
set of this quality in any galaxy, providing an unprece- 
dented opportunity for exploring phase space for forma- 
tional relics. 

We present a basic overview of the GC observations in 
Figure [n^a,b), with a plot of positions in real-space along 
with a phase-space plot of velocity vs. galactocentric ra- 
dius. It is important to note that most of the spatial 
substructure apparent in the GC data is just a prod- 
uct of inhomogeneous spectroscopic sampling: e.g., the 
"gap" at a radius of ~ 7 arcmin, and elongated exten- 
sions to the far East and Northwest. Our analyses will 
be relatively insensitive to such inhomogeneities since we 
will focus on the distribution of velocities as a function of 
galactocentric radius. Furthermore, the sampling effects 
will be built in to our statistical tests for substructure. 

Most of the GCs show a broad spread of velocities re- 
flecting a dynamically relaxed system, but there are two 
sharper features that stand out in phase-space. The first 
is a small elongated "stream" of GCs at 150 kpc, as- 
sociated with the known stellar filament. The second is 
a large chevron-shaped overdensity between 50 and 100 

* Our sample of "GCs" includes some mysterious objects termed 
ultra-compact dwarfs, which have sizes of ~ 10 pc or more, in 
contrast to classical GCs with sizes of a few pc. Whether or not 
there are two or more classes of object, or a continuum of GCs with 
different propert ies, is a major scie ntific question that we cannot 
resolve here (see lBrodie et aLll20111) . Therefore although we know 
for a small subset of our spectroscopic sample that the sizes are 
unusually large, we retain these in any case as tracer particles in 
the halo (excluding only S923 which has a peculiar morphology). 



kpc, which we denote the "shell" . It is an enormous 
structure spread out over ~ 10^ kpc^ on the sky, and is 
not identifiable with any previously known photometric 
feature (Figure[lja)). The statistical significance of these 
two substructures will be demonstrated in Sections |3] and 
m where we also find independent support from their dis- 
tinct GC metallicity distributions. 

Both features are most naturally explained as the 
shredded remnants of infalling galaxies, which we will 
investigate in Section [5] using numerical simulations. To 
orient the discussion of the rest of the paper, we present 
illustrative results from these simulations in Figure [21 
which broadly reproduce some of the general features of 
the data. 

3. INNER SHELL ANALYSIS 

Here we characterize the properties of the large shell 
identified in Figure [Ub), and attempt to establish its 
statistical significance. We exclude the outer-stream data 
set from these analyses, thereby focusing on the main 
part of the data set which has relatively unbiased spatial 
sampling. 

In Section [3. II we compare the use of different conven- 
tions for radial distance. We carry out a group-finding 
analysis and quantify the entropy in the phase-space data 
in Sections l3.2l and l3.3[ and perform maximum-likelihood 
fitting of simple shell models in Section 13.41 In Sec- 
tion [23] we analyze the shell GC color distributions, and 
in Section [3.61 consider the implications for the shell pro- 
genitor. 

3.1. Radius convention 

One of the greatest challenges in detecting and ana- 
lyzing substructures in extragalactic systems is decoding 
the multi-dimensional information contained in observa- 
tions. The dynamical definition of a substructure is that 
it occupies a distinct region in some phase-space of phys- 
ical parameters that are linked to orbital trajectories. 
This could be a space defined by the integrals of motion, 
such as the energy and angular momentum in a spherical 
gravitational potential, or it could be the six-dimensional 
observational phase-space of position and velocity. 

To pick out a substructure from a background of unre- 
lated structures, it helps to link the observations to or- 
bital solutions. This can be tricky enough with 6-D data 
as available for some objects within the Milky Way, but 
in distant galaxies, only three dimensions are normally 
observable (two of position, one of velocity). This "pro- 
jection" of phase space typically smears out the orbital 
tracks and can make substructure inferences degenerate. 

If one already knew a priori the orbital characteristics 
of a substructure, one could apply an optimal coordinate 
transformation to the observations that would maximize 
the contrast in a "reduced" phase-space plot. One of 
the most basic reductions to try is converting the two 
spatial positions to a galactocentric radius, since in the 
limit of a spherical potential and either a radial or a 
circular orbit, the particles will outline a characteristic 
chevron-shaped track in the phase- space of radius and 
line-of-sight-velocity (cf. Figure 1 of lRix et aLlfTool . 

Real galaxies are not perfectly spherical, and in dy- 
namical modeling there is a standard modification to the 
above transformation that corrects for flattening to first- 
order. This is to use an ellipse-based circular-equivalent 
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Fig. 1. Overview of spectroscopic GC observations around M87. Orange circles show confirmed GCs, red crosses show five centrally- 
located low-luminosity elliptical galaxies, and candidate substructure members are highlighted with blue and purple. These are an inner 
shell (at radii of ~ 10-20 arcmin) and an out er st ream (at ~ 35 arcmin), where the most likely members are identified by maximum- 
likelihood fitting to simple models (see Section |3.4| I. (a): Data locations in positional space, overplotted on a deep optical image (down 
to a V-band surface brightness of fiy ~ 28.5 mag arcsec~'^ : IJanowiecki et al.ll20l6l) . with a 100 kpc (21 arcmin) scale illustrated by a bar 
with arrows, (b): The phase space of line-of-sight velocity vs. galactocentric radius, with the dot-dashed horizontal line representing the 
systemic velocity of M87 (1307 km s~^). The radius plotted is the intermediate-axis radius (see Section 13.111 . The velocity uncertainties 
are n ot shown, but are generally smaller than the point sizes in the plot. The satellite ga laxy velocities a re fro m Hypercat IjPaturel et al.l 
\2003 ). where we note that the value for NGC 4486A was recently dramatically revised bv lPrugniel et al.l 1)201 11 ) . 
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Fig. 2. — Simulations of galaxies falling into an idealized central cluster potential, As in Figure [T] both positional space (a) and phase 
space (b) are shown. In each simulation, subsamples of 100 and 10* particles are plotted as large and small dots, respectively. The red 
X symbols show the locations of the progenitor nuclei. Two independent simulations are superimposed, representing initial apocentric 
distances of 90 and 200 kpc (using small black and large orange dots, and small gray and large purple dots, respectively). These two cases 
are intended to be qualitative analogues to the M87 shell and stream (Figure[T]l, with the 0.2% particle subsamples as "GCs" (with small 
measurement errors added). The snapshots shown correspond to 3.5 Gyr (~ 10-20 dynamical times) after the initial apocenters. In the 
first case, there are sharp features in phase space even when the tidal debris is well mixed in positional space (which happens more rapidly 
for a smaller initial apocenter). See Section [5] for further details. 
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Fig. 3. — Radius-velocity phase-space of GCs in M87, omitting the objects in the outer stream region, where different conventions for 
the galactocentric radius are used, (a): simple radius, (b): ellipse-based radius (see Figure (Hb), and the main text for details). The 
ellipse-based radius causes the shell-like feature to appear more distinct. 
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or "intermediate-axis" radius. If x and y are projected 
coordinates along a galaxy's major and minor axes, re- 
spectively, then the simple projected radius is: 



y 



The equivalent radius, given an axis ratio q, is: 



R-n 



' gx^ H . 



(1) 



(2) 



This convention allows one to conserve the area (and 
approximately the mass) enclosed within an ellipti- 
cal isophote when transforming to a circularized one- 
dimensional model. It is also equivalent to recovering 
a constant radius r in three-dimensions for an inclined 
thin disk. 

The next question is what value for q to use, as both 
the flattening of the observable objects and of the unob- 
served gravitational potential could be important. Un- 
fortunately, in the outskirts of M87, neither of the flat- 
tenings is well constrained (i.e., for the GC system or for 
the potential which is dominated by dark matter). It was 
found in S+11 (figure 7) that the GC system has a typi- 
cal q ^ 0.75 at Rm ^ 2'-10' (^ 5-50 kpc), while possibly 
decreasing with radius. At larger radii, the measurement 
is more difficult, but we have used a maximum-likelihood 
technique to estimate a possible range oi q ^ 0.55-0.75, 
at i?,„ - 15'~20' (~ 70-100 kpc). 

Our constraints on q{Rm) for the G Cs are consistent 
with t he profile for the stars derived by iKormendv et all 
(|2009f ). where q declines from ~ 0.9 at Rm ~ 1' to ~ 
0.55 outside Rm ^ 10'. We therefore adopt these stellar 
results for our default model for the GCs, for the sake of 
having smooth, plausible profiles of q and position angle 
with radius 

Figure [3] then compares radius-velocity phase-space 
plots for the M87 GCs, using both the simple and the 
elliptical radius. The cold shell at ~ 50 kpc appears con- 
siderably sharper when using the elliptical radius, and 
we suspect that it may be easier to find substructures in 
other data sets if the Rm convention is used. 

We will adopt this Rm transformation as our default 
for the analyses used in the rest of this paper, while issu- 
ing some caveats. One is that it is not entirely clear from 
a theoretical standpoint why the transformation should 
work: an accretion event generally defines its own orbital 
plane that will not necessarily align with the apparent 
shape of the rest of the galaxy. It could be that we are 
seeing an effect related to the disrupted accretor being 
highly spread out in the host potential. 

The other concern from an observational standpoint is 
that the q{Rm) profile for the GCs is not well determined, 
and fairly small changes in the radius-transformation (of 

^ A more recent study of M87 with deeper surface photome- 
try found that the stellar-Hght flattening becomes even stronger 
at la rge radii, reaching q ~ 0.4 at Rm ~ 20' (.Janowiecki et al., 
120101 ). This is clearly inconsistent with the flattening of the GC 
system, and it would be inappropriate to use this profile for our 
analyses. Both the Kormendy and the Janowiecki profiles agree on 
a position angle of ~ 150° at large radii, while our GC analysis 
is more suggestive of ~ 105°-135°. We do not ascribe much cre- 
dence to this inconsistency, given the difficulties with background 
contamination and spatial incompleteness at large radii, and given 
the much better consistency between stars and GCs at small radii. 
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Fig. 4. — Phase-space group finding in a mock data set con- 
structed from a simulated accretion event (see Figure[2lb)) super- 
imposed on a randomized population of GCs. The genuine shell 
objects are shown as black circles, with velocities that differ slightly 
from the mock data set because of simulated observational errors. 
The general morphology of the substructure is recovered by the 
group-finding algorithm (blue dots). Compare application to the 
real data in Figure nTb). 

^1') would noticeably affect the inferred morphology of 
the substructure, causing it to lose some of its appar- 
ent sharpness in phase-space. Such "blurring" would be 
expected if there is a genuine cold feature, and if one 
deviates from an optimal coordinate transformation as 
discussed above (we have verified this effect using mock 
data simulations that we will discuss later). However, 
we consider it more difficult to artificially sharpen phase- 
space features through a transformation error. 

In summary, the choice of the radius in the phase-space 
diagram is a source of systematic uncertainty that is dif- 
ficult to quantify. Where relevant in our analyses, we 
will also consider the possibility that the substructure is 
more amorphous than in our default model (e.g., as in 
Figure a)). 

3.2. Group finding 

To ascertain the veracity of the shell, rather than rely 
only on the notoriously misled human eye which is good 
at creating patterns out of random noise, we examine sev- 
eral different types of models and statistical tests. This 
is an aspect that has been underdeveloped in previous 
work on GC kinematics around early-type galaxies. 

Our first step is to use a group-finding algorithm, 
based on a standard friends-of-friend s approach (e.g., 
iPerrett et afl [20031: iRlIdick et al.|[2009l ). We define a dis- 
tance between objects of: 



A.S 



' A logio R„ 

WR 



Av 



(3) 



where wu and Wy are radius and velocity scale factors 
that act inversely as weights. Note that our adopted two- 
dimensional phase-space metric differs from the three- 
dimensional {x, y, v) metric commonly used, which is 



8 



Romanowsky et al. 



50 



Rr. [kpc] 

100 150 200 



50 



100 150 



200 



50 



100 150 



200 



2500 



2000 



c/^ 1 500 



^1000 



I I I I I I I I I I I I — I I I I I I I I I I 

real data ■ 

-• % • • • 

-• • • — 

Ay • • , • 

fe^-.-n-i •/ .7. . . , 

••L * « — •- 



500 











I I I I I I I I I I I I I I I I I I I I I I I 



mock data ■ 




9 



•t 



b :: 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

mock data ■ 



m • 



9 • 



5* • •• • 




I I I I I I I I I I I I I I I I I I I I I I I 



10 20 30 40 
R rarcmin] 



10 20 30 40 
R [a rem in] 



10 20 30 40 
R rarcmin] 



Fig. 5. — Group-finding of GCs in phase space. Orange dots show the parent GC population, while green dots show the largest identified 
groups in each data set. (a): the real M87 data (where the outer stream region is again omitted), (b): the mock data set generated from 
an unstructured model, with the largest "group" marked, (c): a typical mock data set. The false group detections typically have many 
fewer objects than in the M87 data, and lack the clear, fairly symmetric shell-like morphology. 



suited for very localized substructures that have had lit- 
tle mixing, but we find that it does not work well for 
structures like shells and long streams that have elon- 
gated "diagonal" morphologies in phase-space. Using the 
3-D metric may have caused large-scale substructures to 
be missed in the past. 

A group is defined as a set of objects separated by less 
than a linking length A in phase space, which is in turn 
defined relative to the mean interparticle spacing (As) by 
a dimensionless parameter A. We then have the criterion: 



|As| < A = A(As) 



(4) 



The next step is to choose the free parameters 
{wji,Wv, X) in the algorithm, where the objective is to 
maximize the sensitivity to real features in the data while 
minimizing the incidence of false groups and linkages. As 
a starting point, we construct mock data sets based on 
simulations of an accretion event (which we will discuss 
in detail later) . To represent the shell, we select 47 parti- 
cles, a number that is chosen to mimic the real data set, 
in the sense that the simulation has 27 particles in the 
same radial region where we find a shell in the real data 
(which has an estimated 27 GCs as we will see later). 

We then fill out the rest of the mock data set by adding 
422 GCs drawn from a Gaussian velocity distribution 
with a = 450 km s~^, and with positions correspond- 
ing to those in the real data set, but randomized by 
20% in order to reduce any lingering spatial substruc- 
tures. Throughout this paper, we will construct similar, 
smooth "background" GC populations, which would be 
unrealistic at some level if the halo were rich in substruc- 
tures of various masses and stages of dissolution. In the 
future, one could consider constructing background dis- 
tributions based on a range of theoretical models for halo 
accretion histories, but for the purposes of this paper, we 
focus primarily on the null hypothesis of a halo with no 



substructure. A caveat to then keep in mind is that the 
large "shell" which we recover could in principle consist 
of two or more distinct substructures that overlap in the 
observable phase-space and conspire to masquerade as a 
single structure. 

Note that the value of tr = 450 km s^^ is larger than 
the ^ 300 km s~^ of the real data, which is done to 
be consistent with the simulations' circular velocity and 
larger shell velocity width. We will later rescale the re- 
sulting group-finding parameters to suit the real data. 
The exact velocities of the mock data set are also con- 
volved with small errors corresponding to the observa- 
tional uncertainties in the real data. 

We find optimum parameters for the simulated data 
set of wr = 0.035 dex, Wy = 50 km s~-^, A = 0.07, with 
acceptable variations around these values at the ^ 20% 
level. The group-finding results are shown in Figure |4j 
a group of 70 objects is identified, which reproduces the 
broad morphology of the shell's limb (cf. Figure Efb)) 
and recovers many of the genuine shell members along 
with many false associations. This exercise provides 
some confidence in the simple group- finding approach for 
recovering cold phase-space substructures that comprise 
only ~ 10% of a sample. However, the individual objects 
assigned to the substructure can include a high degree of 
contamination. This fact, combined with the potential 
sensitivity of the group-finding parameters to the details 
of the data sampling and of the substructure properties, 
motivates exploring alternative algorithms (which we do 
below) as well as the future development of more efficient 
phase-space group-finders. 

Now after rescaling the foregoing simulation-based pa- 
rameters, we adopt the following for use with the real 
data: wr = 0.035 dex, Wy = 30 km s^\ A = 0.07. The 
group-finding algorithm then recovers a shell-like sub- 
structure with 68 members, around half of which may be 
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Fig. 6. — Line-of-sight velocity distribution of GCs in the shell region around M87 (Rm=10.5— 20 arcmin). The systemic velocity of 
1307 km s"'^ has been subtracted from the data. The filled histograms show the data, while the long- and short-dashed curves illustrate 
single and double Gaussian model fits, respectively (with additional smoothing to account for measurement errors). The panels show 
different radial regions, (a): entire region. A mock data set is also shown as a dotted histogram, (b) and (c): finer radial bins as labeled. 
The shell-like substructure is clearly visible as an extra, low-dispersion component in the overall velocity distribution. At smaller radii, it 
appears as a symmetric double peak, and at larger radii as a narrow single peak. 

false associations, which we illustrate with Figure [5{a). 
Given the uncertainty in the radius convention discussed 
in Section rS-H we also carry out the same procedure while 
using the simple radius _Rp. We find in this case a group 
with 80 members, demonstrating that our detection of 
a large phase-space substructure in M87 is not strongly 
sensitive to the radius convention. 

To further evaluate the statistical significance of this 
shell, we carry out Monte Carlo simulations of underlying 
smooth data sets. We generate 100 mock GC data sets, 
using in each case the same radial locations as in the real 
data, but substituting a random velocity drawn from a 
Gaussian distribution with a — 324 km s^^ (which is the 
dispersion found in the real data). 

We then run the group-finding algorithm with the same 
parameters as for the real data. The median largest 
group size found (outside the crowded central regions 
where the group-finding parameters are not optimized) 
has 30 members, and the largest group of all has 61 mem- 
bers. Examples of these groups are shown in panels (b) 
and (c) of Figure [H alongside the real data in panel (a). 
We thus exclude finding a 68-member group by chance 
at the 99% level. This test is somewhat conservative in 
using the radial locations of the real GCs in the mock 
data sets, which may already have some dumpiness. 

If we were to include the data from the outer stream 
region, a stream-like phase-space feature would indeed be 
picked up by the group-finding algorithm, but so would 
be 6 other small "substructures" that (based on the 
Monte Carlo simulations above) are probably not real. 
The outer stream is analyzed separately in a later sec- 
tion, using a slightly different algorithm because of the 
much smaller number counts and the importance of two- 
dimensional spatial clumping. 



where Nj is the number of objects in each partition 
(jHelmi et al.lll999[ ). We consider the 4-20 arcmm region 
of interest and adopt a gridded partition of (radius by 
velocity) with size (2 arcmin x 80 km s^^). We measure 
iS* and compare it to 1000 mock data realizations using 
Gaussian velocity distributions. We find that this region 
is significantly substructured at the 1 a level. This is a 
conservative test because again we have used radii from 
the real data, and because the rectangular partitioning 
is not optimized to pick up extended diagonal features 
as the data appear to show. 

3.4. Maximum likelihood modeling 

We next carry out a set of parametric model fits to the 
data in the shell region, which both allows us to evaluate 
the significance of the substructure and to quantify its 
properties. We will in general focus on the Rm = 10.5- 
20 arcmin (50-95 kpc) region where the shell is most 
apparent by eye. 

The basic idea is that given a model for the line-of- 
sight velocity distribution (LOSVD) dL/dv at a given 
measurement location Ri, we quantify the likelihood C 
of a velocity measurement Vi ± Avf. 



C{vi,Ri) 



dL 
dv 



(v,Ri 



dv 



(6) 



(IRo manowskv fc Kochane3l2001[) . The likelihood can be 
considered as equivalent to a statistic after the oper- 
ation — 21n£. For a Gaussian model dL/dv with an in- 
trinsic dispersion a centered on a systemic velocity Vgys, 
the maximum likelihood model fitting involves minimiz- 
ing the simple function: 



2 In/: = 



usys 



In [cr^ + {Av,f 



(7) 



3.3. Entropy 

The next test uses a very general entropy metric S that 
measures the dumpiness of data in a set of partitions: 



S = - ^N,\nN,. 



(5) 



(|Nolthenius fc Fordl[T986l : iBergond et al.ll2006[ ). 

Before constructing detailed models of the system in 
phase-space, we consider the simplest possible models, 
lumping all of the shell-region velocity data together, 
with no spatial information. Using the likelihood for- 
malism above, we fit a single Gaussian to these data and 



10 



Romanowsky et al. 



2500 



10 20 30 40 

1 — r—TT — |— 1 — I — n — I — I — I — rn — | — r—rn — r-| — i i i 



10 20 30 40 



2000 ;jv .: . . • 

• • • 



,1500 

000 

500 

2500 
2000 
,1500 

i'looo 

500 




-<«>£•.•.• •„ Ts ••• • . 



a 



real '.. tapered Gaussian 



I I I I ' I I I I I I I I I I I I I I I 



C 



* . • • 



• » • * • 
cC*. • •» .* . • 



_L 



_L 



chevron 

J I I I lJ I I I 



b 



^ * . • • • • 

a" • •• • 

• • t.a mm* * 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



d 



* . • . 




wedge 

J I I I I !_] I I I I I I lJ I I I lJ I I I 



10 20 30 40 
[arcmin] 



10 



20 



30 



40 



[arcmin] 



Fig. 7. — Simple models for GC substructure in phase-space, (a): the real data. Vertical dashed lines show the radial range where 
the models are fitted, (b): mock data set using a tapered double-Gaussian model, (c): mock data set using a Gaussian+chevron model, 
(d); mock data set using a Gaussian-!- wedge model. Some single-Gaussian models can also be seen in Figure [5lb,c) for comparison. The 
real data appear to have a phase-space morphology intermediate to these various two-component models, but closest to the chevron-based 
model. 

the substructure, and is insensitive to the adopted radius 
convention fSection 13.11) . since we have not used radius 
information except to define the overah radial bin, and 
the results are also very similar if we define this bin using 
the circular radius Rp. 

Next, given the radial dependence of the shell that is 
visible by eye in Figure E^b), we examine the LOSVD 
in two radial bins: one at the "mouth" of the substruc- 
ture and one at the "apex" (Figure Ufa, b)). We again 
fit double-Gaussian models, allowing for a split around 
Vsys for the cold component. In the first bin, we find a 
dispersion of 23ti| km s'^ with /s = 0.30j:!]:J^ and a 
peak separation of 214 ± 22 km s~^. In the second, we 
find /s = 0.12~^q'q'^, and a dispersion consistent with zero 
with a 68% upper limit of 20 km s^^. 

These binned results together suggest a cold substruc- 
ture with fs ^ 0.2 and a dispersion of ~ 15 km s~^. The 
latter result is similar to our rms velocity measurement 
errors, which can be challenging to determine accurately, 
and consequently a dispersion anywhere in the range 0- 
35 km s^^ is plausible. 

We now model all of the shell-region data jointly in 
phase space. Each model consists of a fraction fs of 



find a dispersion of 315 km s~^. The implied LOSVD for 
this model is shown along with the data in Figure |6ja). 

It is apparent that a single Gaussian is not a good rep- 
resentation of the data, so we next try a double- Gaussian 
model. We find a best fit where the hot component has 
a dispersion of 361 ±34 km s~^, and the cold component 
comprises a fraction fs = 0.25 ± 0.12 of the GCs, with a 
dispersion of 86 ± 34 km s~^ (where Monte Carlo simu- 
lated data sets were used to estimate these uncertainties 
and to correct for slight fitting biases) . We show a mock 
data set of this solution in Figure ID^a) , to illustrate the 
level at which noise can influence the LOSVD. 

By eye, the second Gaussian appears required for ad- 
equately fltting the data, and to quantify this statement 
we carry out a Monte Carlo analysis using a series of 
single- Gaussian mock data sets. We find that fs values 
as high as 0.24 are fitted by chance in only 15 out of 1000 
simulations, while most of these false second components 
are hot, with only 1 out of 1000 having a dispersion as 
low as ^ 80 km s~^. We conclude that there is defi- 
nitely a secondary cold component, with a dispersion of 
120 km s"^ at the most, comprising 13%-37% of the 
GC data. It is also important to note that this con- 
clusion does not depend on the detailed morphology of 
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Fig. 8. — GC line-of-sight velocity distributions, after projection with a diagonal velocity model (see main text), smoothed by 30 km s~^. 
The solid curves shows the overall distribution in each panel, with sample subpopulations shown as dotted curves. The vertical dashed lines 
show lb 160 km s~^ model reference points, (a): real data around M87 from the Rm, = 10.5—20 arcmin region, with the intermediate-color 
subpopulation also shown as a dotted curve, (b): mock data set comprised of a wide Gaussian and a second narrow Gaussian with a tapered 
radial dependence (which is shown as the dotted curve), (c): Gaussian-|-chevron mock data set, with the chevron subsample also shown, 
(d): Gaussian+wedge mock data set, with the wedge subsample also shown. Only the chevron-based data set has a sharp double-peaked 
velocity distribution as in the real data, after projecting with the diagonal model. 



substructure, and a Gaussian component with fraction 
(1 — /s) and dispersion cti. The substructure LOSVD is 
parametrized in one of three ways: as a single Gaussian, a 
nested double Gaussian, or a step function. These mod- 
els are all combined with a linear trend of LOSVD width 
with radius, producing a tapered Gaussian, a chevron, 
or a filled wedge, respectively. The characteristic width 
vs. radius is: 



±v' X (R-R^), 



(8) 



and i?x is the apex radius, 
is set as a fixed parameter 



where v' is the slope 
The systemic velocity 
Ways = 1307 km s~^. 

We fit the various model alternatives to the shell re- 
gion and find fs ^ 0.22, 0.23, and 0.26 for the tapered- 
Gaussian, chevron, and wedge models, respectively. Ex- 
amples of mock data generated from the best model fits 
are shown in Figure [71 By eye, the real data suggest a 
feature that is intermediate between the different mod- 
els, but most closely related to the chevron, owing to the 
strong "edge" seen in the substructure. More quantita- 



tively, Monte Carlo modeling of likelihood fits is unfor- 
tunately not very informative about which model is the 
best representation of the substructure. We find that any 
of the alternative models can readily fit mock data based 
on the other models. 

Instead we turn to a somewhat less parametric rep- 
resentation of the data, using simply a diagonal model 
(with a characteristic slope v' = 17.3 km arcmin"^) 
to project each data point to a velocity at a radius of 
Rm — 10 arcmin. This exercise allows us to collapse all 
of the position- velocity phase-space data into single one- 
dimensional LOSVDs analogous to Figure [51 comparing 
real data to simple models. 

The results are shown in Figure [51 where a 30 km s~^ 
Gaussian smoothing kernel has been applied in order to 
capture the shape of the distribution while reducing Pois- 
son noise. Given an apex in all cases of i?x — 19 ar- 
cmin, the model projects to characteristic velocities of 
'^mod ^ 160 km s~^ (relative to fays) at 10 arcmin. Ap- 
plying the projection to mock data, two peaks are always 
generated but other properties are more model depen- 
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dent: the peak widths, locations relative to Umod, and 
dips in between. In particular, a chevron model when 
added to a Gaussian produces the most distinct, narrow 
peaks close to fmod- 

We find by inspection of Figure H] that the real data 
most closely resemble the chevron-based model, with nar- 
row, well-defined peaks near Vaiod- This plot also indi- 
cates a two-sided rather than one-sided chevron (see also 
Figure |6l). There is furthermore an additional narrow 
feature near —400 km s~^ that appears stronger than a 
random fluctuation and suggests the presence of another 
sharp edge in phase-space terminating at i?„i ~ 30 ar- 
cmin (see also Figures [31(b) and|6l), i.e., perhaps tracing 
a second, fainter shell. 

Although our overall conclusion from LOSVD analyses 
is that a chevron is the best representation for the shell, 
it is difficult to be certain given the available data and 
the uncertainties in the coordinate transformations to 
Rm ■ The real shell probably has a morphology somewhat 
intermediate between the idealized models studied here, 
and also shows hints of being asymmetric with respect 
to Vsys- For now, we report the best- fit parameters for 
the favored chevron model, with uncertainties based on 
Monte Carlo simulations: fs = 0.23 ± 0.06, cri = 19 ± 
10 kms-i, v' = 15.2±1.5 km s'^ arcmin^i, a2 = 358± 
27 km s"^ Of the 119 confirmed GCs in the i?,„ = 50- 
95 kpc region, we thus estimate that 21-35 of them are in 
the shell. We highlight the 27 most likely shell members 
based on this chevron model in Figure [Ijb) . 

Note that the quoted error bars on these fit parameters 
are the statistical uncertainties for an assumed simple 
chevron model. The systematic uncertainties from con- 
sidering alternative models, and from allowing for differ- 
ent radius conventions (Section l3.ip , are difficult to quan- 
tify. Above we reported also a maximal-ignorance model 
that disregards any spatial tapering of the substructure 
by fitting a simple double-Gaussian, and returns a very 
conservative upper limit of 120 km s^^ on the dispersion. 

3.5. Color and metallicity distributions 

One great advantage of using GCs as phase-space trac- 
ers is that they provide additional measurable proper- 
ties beyond position and velocity. Because each GC 
is to a good approximation comprised of a homoge- 
neous stellar population with a single age and metallic- 
ity, the principles of stcUar-subpopulation tagging used 
in the Milky Way (Hclmi 2008) can be extended to more 
distant systems. In particular, GCs are now known 
to be general ly very old object s (ages of ~ 10 Gyr), 
both in M87 (fCohen et al.l 119980 and in other galaxies 
(|Brodie fc Stradedl2006D ." Therefore their optical colors 
are dominated by metallicity variations rather than by 
ages, and we can use color as a good metallicity proxy. 

Rather than solely relying on literature findings on this 
point, we consider the color-metallicity link in our own 
data set. From the DEIMOS spectra, there are 47 objects 
with high enough S/N to estimate iron-based metallici- 
ties [Fe/H] via the CaT line strengths, using techniques 
that we ha ve developed earlie r with the same instrumen- 
tal set-up (|Foster et al.ll2010l 120111 ). Approximate line- 
strength uncertainties are estimated from the spectro- 
scopic signal-to-noise ratio, after calibration from Monte 
Carlo simulated observations of model spectra. 

We provide these metallicity estimates in Tabled and 



in Figure plot them vs. dereddened colors for the 
same GCs where {g — i)o data are available. We find 
a strong correlation between metallicity and colo r that 
agree s roughly with model predictions ( Vazdckis et al.l 
120031 ) . Understanding the scatter in the data will require 
more detailed investigation. 

TABLE 2 

Cat measurements of M87 globular, clusters 



ID R.A. Decl. io {g - i)o CaT 

[J2000] [J2000] [A] 



Note. — Note that the (g — i)o color uncertainties are based 
on shot-noise only, and are useful for relative comparisons of 
colors in this dataset. The absolute uncertainties from calibra- 
tion, reddening correction, background estimation, etc., will be 
larger. The equivalent iron metallicities shown on the right-hand 
axis of Figure [9l a) are derived using the expression: [Fe/H] ~ 
0.438xCaT-3.641. 

This test does not prove unequivocally that these col- 
ors are a perfect proxy for metallicity, given various unre- 
solved issues about CaT-based metallicities, and the pos- 
sibility of some residual age effects on the colors. How- 
ever, it does provide assurance that the current conven- 
tional wisdom about GC color as a reliable metallicity 



S731 

S878 

S952 

S1007 

S1074 

S1265 

H20493 

H21863 

H25785 

H26690 

H27496 

H27916 

H28411 

H28415 

H28866 

H30757 

H30772 

H31134 

T13456 

T13500 

T13559 

T13569 

T13571 

T13589 

T13593 

T13609 

T13611 

T13623 

T13642 

T13648 

T13899 

T13901 

T14060 

T14108 

T14133 

T14149 

T14228 

T15685 

T15777 

T15863 

T15886 

T15900 

T15908 

T16080 

T17211 

VUCD6 

VUCD8 



187.72452 
187.70708 
187.69952 
187.69344 
187.68856 
187.67009 
187.87314 
187.87253 
187.86159 
187.72956 
187.80753 
187.71521 
187.70221 
187.78343 
187.75563 
187.71996 
187.74191 
187.72735 
187.92093 
187.93052 
187.93331 
187.98853 
187.89063 
187.92350 
187.90345 
187.96307 
187.92010 
187.90570 
187.94181 
187.98298 
187.84033 
187.87673 
187.87089 
187.83846 
187.84645 
187.87374 
187.84434 
188.04605 
188.14226 
188.17885 
188.15205 
188.21484 
188.15419 
188.27954 
188.33727 
187.86816 
188.01813 



12.28682 
12.28149 
12.27727 
12.28971 
12.27776 
12.28303 
12.15391 
12.17069 
12.21268 
12.22270 
12.23123 
12.23610 
12.24179 
12.24165 
12.24666 
12.26691 
12.26728 
12.27091 
12.46766 
12.42406 
12.37800 
12.36469 
12.36288 
12.35406 
12.35247 
12.34319 
12.34146 
12.33556 
12.32521 
12.32301 
12.42620 
12.42616 
12.38884 
12.37504 
12.37052 
12.36463 
12.35020 
12.36726 
12.45638 
12.36688 
12.34920 
12.33844 
12.33126 
12.34479 
12.48077 
12.41766 
12.34176 



20.37 

20.62 

20.58 

19.41 

20.29 

19.82 

20.13 

20.54 

19.88 

19.70 

19.56 

20.45 

20.10 

19.71 

19.96 

20.20 

19.84 

20.04 

21.14 

20.55 

21.87 

21.93 

21.42 

21.71 

21.81 

21.95 

22.44 

21.49 

20.39 

21.79 

22.38 

21.83 

20.83 

19.82 

20.69 

21.84 

20.44 

22.08 

21.93 

21.36 

22.33 

20.6 

21.19 

21.0 

21.0 

18.63 

18.96 



0.858 
0.816 
0.779 
1.023 
1.120 
0.772 
0.785 
0.977 
0.914 
0.990 
0.801 
0.756 
0.799 
0.925 
0.994 
0.809 
0.906 
0.821 
0.757 
0.767 
0.904 
0.695 
0.801 
0.773 
0.778 
0.729 
0.793 
0.764 
0.869 
0.731 
0.964 
0.853 
0.722 
0.908 
0.980 
0.773 
0.785 
0.730 
0.796 
0.753 
0.743 



± 0.004 
± 0.004 
± 0.004 
± 0.002 
± 0.004 
± 0.003 
± 0.005 
± 0.006 
± 0.002 
± 0.002 
± 0.002 
± 0.004 
± 0.002 
± 0.002 
± 0.004 
± 0.002 
± 0.033 
± 0.003 
± 0.008 
± 0.004 
±0.012 
±0.015 
± 0.008 
±0.012 
±0.011 
±0.014 
±0.022 
±0.011 
± 0.004 
±0.020 
± 0.009 
±0.013 
± 0.007 
± 0.002 
± 0.005 
±0.013 
± 0.004 
±0.017 
±0.013 
± 0.008 
±0.021 



0.781 ±0.007 



0.896 ±0.001 
0.855 ±0.001 



7.05 ±0.31 
7.89 ±0.53 
4.78 ± 0.48 

6.94 ±0.13 
5.89 ±0.30 

5.23 ±0.20 

7.24 ± 0.42 

8.25 ±0.42 

1.95 ±0.22 
7.56 ±0.16 

5.31 ±0.13 

5.46 ±0.37 
5.64 ±0.26 

7.80 ±0.16 

6.13 ±0.20 
6.04 ±0.28 
3.61 ±0.31 

8.22 ±0.23 

5.06 ±0.15 
3.49 ± 0.26 

6.47 ±0.33 
3.77 ±0.25 

5.81 ±0.14 

4.14 ±0.22 
6.12 ±0.26 
6.58 ±0.23 
7.00 ± 0.45 

4.96 ±0.14 

6.32 ±0.02 
3.99 ±0.21 
7.53 ±0.31 
5.56 ±0.20 
4.39 ± 0.06 

6.23 ±0.01 
6.28 ±0.06 
6.42 ±0.31 
5.93 ±0.03 

5.37 ±0.31 
5.41 ±0.19 
5.64 ±0.11 
4.99 ±0.29 

4.20 ± 0.06 

4.38 ±0.08 

6.21 ±0.07 
5.08 ± 0.08 
5.69 ±0.01 
6.06 ±0.01 
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Fig. 9. — Color vs. motallicity for M87 GCs. (a): results from CaT analysis of the subset of our DEIMOS spectra with good enough 
signal-to-noise {S/N > 12) and with {g — i)o data available (44 objects). The left axis shows the direct CaT line strengths, where the error 
bars indicate 68% uncertainties. The right axis shows the equivalent iron metallicity. The color uncertainties in some cases are smaller 
than the point sizes. One additional outlier is off the scale at {g — i)o = 0.91, CaT=1.9 A. (b): color versus metallicity for 148 GCs from 
earlier spectroscopic analysis (Cohen ct al. 1998). Note that while in this case the line-index metallicity conversions below [M/H] ~ —1.8 
may well be non-linear, and above [M/H] ~ the metallicities are not well calibrated, the important issue here is not the exact metallicities 
but rather the monotonicity of the color-metallicity relation, which is demonstrated by both data sets. 



proxy is probably sound in the halo of M87 (i.e., at radii 
of ~ 45-220 kpc). 

Previous spectroscopic work using bluer wavelengths 
at smaller radii 3-40 kpc) also found M87 GC metal- 
licity indices to correlate strongly with U — R color, 
and for the ages t o be almost all larger than 10 Gyr 
(jCohen et al.lll998[ ) . We revisit these data by comparing 
our {g — i)o colors to these authors' metallicity estimates 
in Figure injb), finding a tight correlation between the 
two quantities. Again, this supports the color-metallicity 
link. 

Despite these arguments that color is a good proxy for 
metallicity in our data set, it is important to note that 
our analysis below is independent of this link. Whatever 
the root cause of variations in color, we can still use it as 
an empirical tagging factor. Then for our current pur- 
poses the issue is more a matter of semantics: whether 
we are analyzing chemo- dynamics or chromo-dynamics. 

Our first basic goal is to use color as a completely inde- 
pendent dimension in the data, to see if it supports the 
differences among GCs that emerged from the position- 
velocity phase-space analyses. We plot velocity versus 
color in the shell region in Figure [TUfa). In M87, as 
in galaxies in general, the population of metal-poor GCs 
has a more extended distribution from their host galaxy's 
center than the more metal-rich GCs. Given that all of 
these objects reside in the same gravitational potential, 
and assuming that they have similar orbital anisotropics, 
their spatial differences are expected to be related to their 
velocity dispersions. The metal-rich GCs should gener- 
ally have lower dispersions, and Figure rTOT a') is expected 
to show a systematic decrease in dispersion with color, 
with perhaps a strong transition between the two main 
subpopulations, at (g — i)o ^ 0.9. 



Figure llOf a) does show a kinematic transition with 
color, but not in the way just described. There is an 
overall pattern of a cold component near Vgys embedded 
in a hotter component (cf. Figure [5]), but the cold com- 
ponent does not appear bluewards of {g — i)o ~ 0.77. 
This color is in the middle of the main blue peak for 
M87's GC system, and implies that there is kinemati- 
cal substructure within this peakl3 We have looked at 
the same information using the CaT-based metallicities, 
with consistent results (colder material at CaT ^ 5.3 A 
or [Fe/H] ^ — 1.3 dex), but the statistics are too poor 
to be conclusive (17 GCs in the shell region). 

We next take a different approach, of separating GCs 
into probable shell and non-shell objects, based on 
the phase-space maximum likelihood models constructed 
previously. We then examine the color distributions of 
the two subsamples. For each of the models (tapered 
Gaussian, chevron, wedge), the two color distributions 
are different, with the chevron case shown in Figure fTUf b) 
as an exampleQ Here the shell appears to contain an ex- 
cess population of intermediate-color GCs, peaking at 
(5-z)o~0.90. 

Motivated by this result, we also show the velocity dis- 
tribution of all the intermediate color objects in the shell 
region, with {g — i)o — 0.85-0.95, in Figure HJa). These 
objects do show a double-peaked profile as expected in 
the chevron model, except for an additional low-velocity 
peak that accounts for much of the candidate outer shell. 
We have also fitted the chevron model to the phase-space 
data for the blue, intermediate-color, and red GCs sep- 

This is again a test that is insensitive to the radius transfor- 
mation Rm- 

^ Among our probable shell sample, there are two ultra-compact 
dwarf candidates: VUCDl and H59823. 
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Fig. 10. — Correlations between phase space and color properties of M87 GCs in the Rm = 10.5—20 arcmin shell region, (a): line-of-sight 
velocity versus color. A "cold" velocity component appears at colors redder than {g — i)o ~ 0-77. GCs assigned to the chevron are marked as 
filled orange circles, (b): the chevron and non-chevron subsamples (with 27 and 92 objects, respectively). Histograms show the raw counts 
in color bins of 0.04-mag width. Alternative, unbinned color distributions are shown by the dotted curves, which are produced by smoothing 
the discrete data with a 0.036-mag optimum Gaussian kernel (a standard technique in studies of GC color distributions; Peng ct al. 2006). 
The amplitude of the smoothed chevron distribution has been rescaled to the same integrated number as the non-chevron. The chevron 
GCs appear to have a slightly redder "blue peak" than the non-chevron GCs, as well as a strong intermediate-color peak. 

arately. The resulting parameters for each of these fits 
are mutually consistent, and color does not discriminate 
sharply between shell and non-shell objects, but we do 
find that the shell fraction is highest in the intermediate- 
color population, with fg ^ 0.5. 

We evaluate the statistical significance of the differ- 
ences between the chevron and non-chevron color dis- 
tributions using a Kolmogorov-Smirnov test. We find a 
D value of 0.25, with a 13% probability of occurrence 
by chance. According to the test, the largest difference 
between the distributions comes at the far blue end, sup- 
porting the transition seen in Figure llOf a) . The same is 
true for all three shell models, with the far-blue differ- 
ence being strongest for the tapered-Gaussian case, and 
the intermediate-color peak for the chevron case. The 
similarity of the three models also suggests that these 
conclusions are insensitive to the choice for radius Rm, 
as is also explicitly the case for the test in Figure rTOT a) . 

We have found that metallicity information (primar- 
ily using color as a proxy) provides independent support 
for the presence of the M87 shell. The colors do not yet 
differentiate clearly between the alternative phase-space 
models for the shell, but even our unprecedented data 
set only scratches the surface of what can now be accom- 
plished. It is technically feasible to increase the number 
of measured GC velocities by an order of magnitude, e.g., 
to ~ 1500 objects in the shell region alone. 



3.6. Progenitor inferences from GCs 

Having established the existence of the shell, and some 
of its properties, we now attempt to identify the nature 
of its progenitor, assuming that the accretion and disrup- 
tion of a single galaxy was responsible. Unfortunately, 
we will see that it is difficult to derive a self-consistent 



solution from the basic shell parameters, so we will carry 
out an inventory of different types of constraints to see 
if an overall solution emerges (with an estimate of the 
progenitor's 1^-band stellar luminosity, Ly, as a basic 
goal). 

Since in the previous Section we have just discussed 
the color distribution of the shell GCs, it is natural to try 
using these colors to constrain the progenitor Ly- Em- 
pirica lly, GC colors correlate with host galaxy luminosity 
(e.g., iLarsen et al.l feOOl: Stradcr ct al. 2004; Lotz ct at] 
I2004D . which is conceptually related to the correlations 
between sub structure metallicities and their progenitor 
masses fe.g.. iJohnston et al.ll2008l) . This framework was 
used to estimate the relative masses of acc r eted g alaxies 
in the halo of a nearby SO in lArnold et al.l (|2011 ). 

In the case of M87, we can make use of the correlations 
found from the "ACSVCS" large survey of elliptical and 
lenticul ar galaxies, both giant and dwarf, in the Virgo 
cluster (jPeng et al.ll2006[ ). Ahhough the GC color dis- 
tributions in these galaxies are in general complex, they 
show a general tendency for bimodality, with distinct blue 
and red peaks. The color-luminosity relation is therefore 
reported separately for the blue and red GCs. The un- 
derlying basis for these correlations is thought to be a 
mass-metallicity relation, but again, our analysis does 
not rely on this conclusion. 

We have modeled the color distributions of the chevron 
and non-chevron GCs (Figure [TOT b)) using Kaye's Mix- 
ture Models based on heterosceda stic Gaussian color sub- 
populations (jAshman et al.llT994[ ) . We find marginal sig- 
nificance {p = 0.325) for bimodality among the chevron 
GCs, and high significance (p < 0.001) for the non- 
chevron GCs. Beginning with the non-chevron GCs, we 
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Fig. 11. — Parameter-space of candidate progenitor galaxies for the M87 shell, plotting V-band luminosity, central stellar velocity 
dispersion, and GC numbers. The bands show constraints from modeling the shell GCs, with the narrower purple bands stemming from 
our default chevron-based analysis, and the wider blue bands showing a model that does not incorporate spatial information (yielding 
a conservative upper li mit on the dispers ion). The points with error bars show early- type galaxies in Virgo fr om the A CSVCS survey, 
with A^Qc and Ly from lPeng et al.l I I200SI1 . The dispersions are mostly from HyperCat with some additions from lChiling arian (2009) and 
ri'oloba et al.l 1 120111 ). Points with red X symbols superimposed show four low-luminosity ellipticals near to M87, whose positions in these 
diagrams are suggestive of histories of stripping — as illustrated schematically by the gray arrows. The constraints from the chevron model 
do not match up well with any realistic candidate progenitor galaxies, while the conservative model is more viable (see text for details). 



find that they have a primary peak at (g — i)o ^ 0.77 
and a possible weak secondary peak at — i)o ^ 0.94, 
which correspond to — 2;)o ^ 0.92 and ^ 1.17, or 
[Fe/H] ^ —1.5 and —0.8, respectively (where we have 
made empirical color conversions using a large set of M87 
GCs that overlap between ACS and CFHT imaging). 

We next assume for simplicity that the GC color peaks 
in the surviving Virgo early-type galaxies are represen- 
tative of the disrupted progenitors that built up the M87 
haloQ Considering first the peak locations of the non- 
chevron GCs and comparing them to the non -parametric 
peak s determined for the ACSVCS galaxies (jPeng et alJ 
l200l) . we find the former to be fairly atypical when con- 
sidered jointly (a puzzle that would be worsened if the 
chevron GCs were included). This makes the progenitor 
inference unclear, perhaps owing to the superposition of 
color peaks from multiple progenitors, but our best guess 
would be a typical progenitor magnitude of My ~ — 18 

[log(Ly/Ly,0) ^9.5]. 

Approximating the T^-band stellar mass-to-light ratios 
as roughly constant, the implication would be that the 
accretion events which built up the bulk of the M87 halo 

® There is evidence that galaxies accreted during earlier 
epochs had systematic ally redder GC colors at a fixed mass 
flBrodie &: Stradeill2006l ). which would mean that our luminosities 
for the typical M87 halo progenitors will be overestimated. 



had characteristic stellar- mass merger ratios of ~ 1:30. 
This may be consistent with current ideas of massive 
galaxies' stellar e nvelopes being a ssembled through mi- 
nor mergers (e.g.. lOser eFall [20111) . 

Turning now to the chevron GCs, they have a pri- 
mary peak at {g — i)o 0.79 and a secondary peak at 
{g — i)o ^ 0.90. corresponding to {g — z)o ^ 0.95 and 

~ 1.11, [Fe/H] 1.4 and -1.0. This red peak location 

is very atypical for the ACSVCS sample, and may have 
been inaccurately estimated through our procedure of 
statistically disentangling the chevron and non-chevron 
GCs. We therefore consider the color constraint on the 
shell progenitor luminosity to be indeterminate. 

The other two major clues about the shell origin are 
its number of GCs, NqCi E^nd their mutual velocity dis- 
persion, iJgc- For both of these properties, we use the 
results of maximum likelihood models from Section 13.41 
While our default constraint comes from the "chevron" 
model, we will also account for the systematic uncer- 
tainty in the shell morphology (e.g., owing to the radius 
conventions discussed in Section 13. ip by considering an 
alternative scenario of maximal ignorance of the mor- 
phology. For this case we adopt the double- Gaussian 
analysis of Sect ion [3T4l which leads to a much larger value 
of CTGC (50-120 km s^^ rather than 10-30 km s'^). 

To make use of the A^gc and ctgc constraints, we must 
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correct both observational measurements to the relevant 
quantities for the progenitor galaxies. In the case of 
A^GC: we extrapolate the fraction of shell GCs in our 
spectroscopic sample to the much larger sample of unob- 
served GCs in the same region. We have constructed a 
smooth model surface density profile of the M87 GC sys- 
tem, based on both the CFHT photometric data as well 
as spectroscopic data to confirm the non-GC contamina- 
tion rate (see S-l-11). This model is initially based on 
our photometric limit of i = 22.5 {z ~ 22.4), but to get 
total GC numbers we need to extrapolate to fainter lumi- 
nosities, using an evo lved Schechter funct ion fitted to the 
central GCs of M87 (| Jordan et all 120071 ). This provides 
a scaling factor of 2.760 

Now integrating the revised surface density between 10 
and 20 arcmin where the shell signal is clear, we find a 
total of 3500 ± 400 GCs. Multiplying by the inferred /s 
values for both the chevron and double-Gaussian mod- 
els (see Section 13. 4|) , we estimate A^gc = 800 ^270 and 
870j;|go total GCs in the shell, respectively (integrating 
over all magnitudes). These numbers could be lower if 
we have detected a disproportionate fraction of the shell 
objects by spatial bias in our spectroscopic selection that 
coincides with a spatial clumping of the shell; or higher 
if the shell extends farther in radius and velocity than 
in our simple model (see, e.g.. Figure [2jb)). It is diffi- 
cult to quantify such uncertainties, but we suspect that 
they could involve an additional factor of 2, which 
would accommodate a range of up to ^ 200-2900 GCs. 
If Ngc ~ 800, then the shell accounts for - 20% of the 
GCs around M87 within the 50-100 kpc range, and ~ 5% 
of the total GC system^ 

What kinds of galaxies are observed so many GCs? To 
provide a visual overview of this constraint, we plot the 
number of GCs versus host luminosity for the ACSVCS 
sample in the upper left panel of Figure 1111 finding that 
~ L* luminosities of of log (Ly/Ly^©) ^ 10.1-10.8 are 
implied. 

Next we consider the shell's internal velocity disper- 
sion, which as a first approximation, we assume to be 
similar to the internal velocity dispersion of its progeni- 
tor, as might be expected under the simplest tidal strip- 
ping models. Here the comparison to the ACSVCS galax- 
ies is less straightforward, since very few of these have di- 
rect measurements available of their GC velocity disper- 
sions. Instead, we will adopt a conversion between ctgc 
and the widely-available central stellar velocity disper- 
sion (Tc- Based on the forthcoming kinematical analysis of 
a handful of low- luminosity E/SOs (Pota et al., in prepa- 
ration), we find that the ratio (Tc/cgc varies from galaxy 
to galaxy over a range ^ 1.3-2.3, with a typical value 
of ~ 1.8. We therefore infer CTc 35 ± 25 km s~^ and 
^ 150 ± 100 km s~^ based on the chevron and double- 
Gaussian models, respectively. 

^ With a Gaussian luminosity function, this factor would be 2.50, 
but for direct comparison with A^gc values from the ACSVCS, we 
adopt 2.76. 

To derive the latter value, we have integrated the GC spatial 
density profile out to a radius of 45 arcmin (220 kpc), which is 
where the planetary nebula (PN) kinematics suggest a transition 
between objects bound and not bound to M87 IIDohertv et al.l2009l : 
the kinematics data for the GCs themselves are less extensive at 
large radii and do not yet show any sign of such a transition) . The 
total GC population is then 16200. 



We plot this constraint again versus the ACSVCS sam- 
ple in the bottom left panel of Figure[TTJ Here we see that 
the chevron-based velocity dispersion suggests a dwarf 
E/SO with log(Lv/Lyo) < 9.5EI while the conserva- 
tive upper-limit dispersion allows for luminosities of up 

to log(Ly/Lyo) ~ 10-6. 

For the chevron model, the dispersion constraints on 
luminosity are thus inconsistent with the constraints 
from GC numbers. To make this comparison more di- 
rectly, we show dispersion versus GC numbers in the 
upper right panel of Figure [TlJ The intersections of the 
vertical and horizontal bands mark out possible solutions 
for the progenitor based on the shell analyses — but these 
do not overlap with the properties of known Virgo galax- 
ies. It does not seem plausible that either a dwarf galaxy 
hosted hundreds of GCs, or that an ^ L* galaxy left ex- 
tensive debris with a velocity dispersion of <, 30 km s"^ . 

The conservative model fares much better in this 
respect, as several ACSVCS galaxies are consistent 
with both its Nqq, and ctgc constraints. Some 
examples include NGC 4435 and NGC 4473, with 
\og{Lv / Lv,q) ^ 10.1-10.2. This model seems promis- 
ing, but one should remember that is intended to rep- 
resenting a bracketing case of maximal ignorance about 
the shell morphology, while there does seem to be some 
degree of substructure tapering with radius in phase- 
space, independently of the radius convention. Therefore 
it is worthwhile reviewing any additional lines of evidence 
about the accretion progenitor. 

If the progenitor were very massive, then it should 
have left other clear signatures of its passing, besides 
the phase-space shell. One might expect to see obvi- 
ous disturbances and asymmetries in the light distribu- 
tion of M87, as well as strong merger-induced halo rota- 
tion, which is not observed (S-l-11). If we assume that 
this rules out a mass ratio of ~ 1:3 or higher, it im- 
plies log (Ly/Lyo) < 10.4 for the progenitor, or perhaps 
somewhat higher if dark matter is taken into account. 

Besides the impact of the accretion event on M87, 
the incoming galaxy should have deposited not only 
GCs but field stars, which would leave a large-scale sur- 
face brightness fluctuation across the face of M87. In 
the shell region where the typical surface brightness is 
^.v ~ 27 mag arcsec"^, we can rule out obvious fluctu- 
atio ns at the ^ 0.5 mag ar csec"^ level (see Figure [TJa) 
and iJanowiecki et al.l 120101 ) . We can therefore rule out 
\og{Lv / Lv,q) ^ 9.8 of accreted stars that are dis- 
tributed uniformly over the shell regionQ 

An alternative mechanism for producing a low velocity dis- 
persion is through tidal tails stripp ed out from the cold disk 
of an infalling spiral galaxy (e.g., IHernguist fc Spergell 119921 : 
IHibbard fc Milio£lfT995l) . However, we do not know of any galaxy 
with hundreds of GCs contained in a cold disk. Alternatively, some 
GCs could have been formed from cold gas in such a merger, and 
then deposited on closely associated orbits. Because of the age- 
metallicity degeneracy, the shell GCs might in principle be young 
and metal-rich rather than old and metal-poor. Full analysis of the 
GC ages is beyond the scope of this paper, but our initial inspec- 
tion of the Hectospec data indicates there are no objects as young 
as ~ 1 Gyr. 

This constraint is relevant not only to our default single- 
progenitor scenario, but also to a scenario involving multiple small 
galaxies that together deposited a large number of GCs with cold 
kinematics: e.g., from ~ 5 dwarf E/SOs with log {Ly / Ly q) ~ 9.3 
each. Besides the low probability of accreting such galaxies on 
near-identical trajectories, the associated stellar debris could vio- 
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This luminosity Umit is lower if the stars are confined 
to only one side of M87, as Figure [TJa) suggests; and 
higher if the stars have become spatially decoupled from 
the GCs and reside at smaller radii where they are either 
hard to detect, or are still intact as one of the nearby low- 
luminosity elliptical galaxies. Such spatial segregation is 
plausible at some level since the GCs would be expected 
to inhabit preferentially the outer parts of the incoming 
galaxy, and thus be stripped off before the stars. Note 
that a fluctuation in the surface density of the GCs them- 
selves would be challenging to detect at the ^ 50% level. 

As mentioned previously, there are no major signa- 
tures of interactions previously known for M87. Some 
candidates for large-scale, low-sur face brightness stellar 
features have been id entified (lArp fc Bertolal 119711: 
Weil. Bland-Hawthorn, fc Mahd Il997t IRudick et all 



201(f). but these are at spatially disparate locations 
from the GC shell, and in some cases are likely to be 
manifestations of foreground cirrus. 

More intriguingly, as mentioned above, there are five 
low-luminosity ellipticals within a 30-80 kpc projected 
radius of M87: NGC 4476, NGC 4478, NGC 4486A, 
NGC 4486B, and IC 3443 (see Figures [Hand [HI). Some of 
these have previously been proposed as the stripped rem- 
nants of larger galaxies, based on disturbances in their 
isophotes or unusual stellar colors for their luminosities. 
If any of these galaxies is being actively disrupted, some 
stars and GCs could have already been stripped out and 
deposited within the halo of M87. 

The first four of these galaxies are in the AC SVC S sam- 
ple, so we consider this scenario in more concrete terms 
by marking them in Figure [TT] with red x symbols. Here 
we see from the top two panels that all four galaxies ap- 
pear depleted in GCs relative to the average trend, which 
suggests stripping of their outer regions. The offsets are 
strongest for NGC 4486 A and NGC 4486B, which are 
statistically consi stent with hostin g zero GCs (see fur- 
ther discussion in lPeng et aLll2008[ ). NGC 4486B is also 
a very unusual outlier in the ac-Ly panel, and p ossibly 
in su permassive black hole correlations (Giiltckin et al.l 
12009'), implying that not only its GCs but also 90% of 
its stars have been stripped away. 

The CTc values of NGC 4486A and NGC 4486B sug- 
gest that they originally hosted ~ 70-200 and 100-500 
GCs, respectively. Their phase-space positions in Fig- 
ure [Ijb), where their velocities are ~ 750 km s~^ and 
~ 1550 km s~^, respectively, indicate that NGC 4486B 
could be associated with the shell. NGC 4478 (at 
~ 1400 km s~^) is also consistent with the shell in phase- 
space; it does not appear to be as severely stripped, but 
still could have lost up to ~ 250 GCs. 

The colors of the NGC 4478 GCs also provide one of 
the closer matches to the shell GC colors (this test is not 
possible for NGC 4486A and NGC 4486B since they have 
virtually no GCs remaining). For both NGC 4478 and 
NGC 4486B, the dynamical friction timescales at their 
projected radii are 1-3 Gyr or shorter, reinforcing the 
likelihood that these are objects in the final stages of 
orbital decay toward the center of M87. For both galax- 
ies, there are also outer stellar velocity dispersion mea- 
surements available, making the extrapolation to ctgc 
less uncertain. For NGC 4478 this is 50-100 km s'^ 

late the surface brightnesses constraint. 




Fig. 12. — Closer view of the M87 V-band image from Fig- 
ure [JI a), showing a 23.5' x40' (110x190 kpc) region. Some small 
galaxies in the field are labeled. The shading is logarithmically 
spaced by surface brightness, and arbitrarily colored in order to 
enhance any isophote twisting (e.g., in the case of NGC 4478). 



()Hallidav et aLlllQOll), and for NGC 4486B, it is ^ 100- 
120 kms-i ()Kormendv et al.llT997HSpolaor et al.i l2010). 

The present-day luminosity of NGC 4478, 
\og{Lv/Lv,0) = 9.8, suggests that it may have 
experienced very little stellar stripping (Figure lll[ 
lower- left), which would explain the lack of surface 
brightness variations around M87. For NGC 4486B, the 
implied luminosity at infall of \og{Lv/Lv.Q) ~ 10.1 
would require that much of the stellar debris is hidden 
at small galactocentric radii, which is plausible given 
that this stripped galaxy is currently at a radius of 

30 kpc (note that the missing debris is an issue even 
without considering the GC shell). 

Both of these galaxies appear to be possible shell pro- 
genitors, as their properties are consistent with all of the 
constraints discussed above, within the uncertainties, ex- 
cept for the low velocity dispersion in the chevron model 
(this dispersion could then be considered either as an 
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artifact of the adopted model, or as a product of strip- 
ping dynamics that are not yet understood). In fact, the 
problem could be turned around: given that there are 
stripped galaxies around M87 which must have lost hun- 
dreds of accompanying GCs, should we not expect to see 
some evidence for them in phase-space? 

The multiplicity of candidate shell progenitors is also 
a reminder that our initial single-progenitor assumption 
may be oversimplified. It is indeed possible that the 
phase-space of GCs in the M87 halo contains multiple, 
overlapping, prominent shells and streams More GC 
spectra are clearly needed to better determine the sub- 
structure properties, and to make clearer progenitor in- 
ferences. 

As a final plausibility check, we examine whether or 
not our preferred solution, with an accreted luminos- 
ity of \og{Lv/Lv,Q) 10, (i.e., ~ 0.5 L*, or a - 1:10 
merger), represents a likely event in the context of other 
observational constraints on merger rates in BCGs. As 
discussed in Section [U the details of BCG assembly are 
not well known, which motivates our phase-space study, 
but there are still some plausible upper limits on merger 
rates. For example, M87 cannot have accreted a galaxy 
with \og[Lv / Lv,q) ~ 10 more often than once every 
~ 1.5 Gyr, averaged over a Hubble time, or else the total 
luminosity would be higher than is observed today. 

The calculation is more complicated than this because 
the merger history will be divided over a spectrum of 
mass ratios, and because the merger rate may change 
with time. We therefore turn t o the low-z BCG survey 
carried out by iLiu et all (|2009D , who found that major 
mergers (down to a ~ 1:4 mass ratio) were visible in 1 
case out of every 30. If we assume for simplicity that 
mergers are detectable in phase-space for roughly twice 
as long as in real-space (which is the only place that 
merger timescales will enter into our calculations), then 
recent major merger signatures could be found in ~ 7% 
of BCGs. 

Next, we should like to extrapolate to lower-mass ra- 
tios. The simplest assumption we can make is that the 
accreted objects follow a standard Schechter luminosity 
function, where 1:4 events correspond to ^ L* in the 
case of M87. We then calculate that events down to 1:10 
should be present in ^ 25% of BCGs. Thus it does not 
seem improbable that the shell in M87 traces a ~ 1:10 
merger event, while more massive progenitors would be- 
gin to stretch the limits of plausibility. 

4. OUTER STREAM ANALYSIS 

We next analyze the GCs in the northwest area, around 
"stream A" . Figure [13] summarizes the positional space 
and phase space properties of this area. The known 
stellar stream has a y-band luminosity of (2.4 ± 0.3) x 
lOfLy0 over a ~ 15 x 100 kpc region (jJanowiecki et all 
[2OT0h .' 

It is important to keep in mind that the spectroscopic 
observations were tiled closely around the stellar stream, 
and much of the apparent narrowness in real space of the 
overall GC distribution is an artifact of this bias. The 
fractional return of bona fide GCs from the spectroscopic 
observations was indeed higher "on-stream" than "ofF- 

There may be another, related feature in M87: S-l-11 identified 
a pecuhar velocity ofFset in the metal-rich GCs inside ~ 10 kpc. 



stream" , but the statistics here are relatively poor and a 
more in-depth analysis will be required to consider this 
question further. 

We focus instead on a possible cold phase-space sub- 
structure within the spectroscopic sample, where the hot- 
ter background may be the stream or the general diffuse 
outer GC system of M87. This tentative substructure 
consists of 6 GCs (out of 19 in the sample apart from 
one associated with a nearby galaxy) that are clustered in 
real-space and follow a cold diagonal path in phase space 
(drawn schematically in Figure ITST b')'). Two of these ob- 
jects are a very close pair, separated spatially by only 
2.6 arcsec (0.2 kpc) and in velocity by only 13 km s~^, 
which is consistent with zero within the uncertainties. 

To estimate the chances of this substructure being a 
chance clumping, we pursue friends-of-friends group find- 
ing as previously carried out for the shell (Section 13. 2|) . 
but reverting to a traditional three-dimensional (x, ?/, v) 
phase-space search. In order to pick up the apparent 
group in the real data, we adopt Wx — Wy = 50 arcsec 
(4.0 kpc), w„ = 30 km s~\ A = 0.29. We then simu- 
late mock data sets using the same positions while draw- 
ing velocities from a Gaussian LOSVD with dispersion 
271 km s~^, and find groups with 6 or more members by 
chance only 9% of the time. Again, this analysis is con- 
servative in that the mock data sets start from spatial 
positions that may not be random, but not conserva- 
tive in that the group-finding parameters are optimized 
a posteriori. 

We next carry out maximum-likelihood model fits to 
the data, using a Gaussian plus one-sided chevron model. 
We estimate that A'^^ of the GCs in this region belong to 
the stream, which has a dispersion 

of 33^25 kms-i. This 
model is intended as a characterization of the stream 
properties rather than as a robust test of its reality, since 
the fits are somewhat unstable to the starting conditions. 
The model does not incorporate the full spatial informa- 
tion about the region, and given that the stream objects 
are also clumped in the azimuthal direction, we consider 
it likely that the 5 most closely associated objects seen in 
Figure [l3] comprise the stream. Extrapolating to fainter 
magnitudes from our spectroscopic limit of «o ~ 22.5, we 
estimate a total of ^ 15 GCs in the stream clump. How- 
ever, our spectroscopic survey is incomplete, and there 
might easily be ~ 20-30 stream GCs. 

As with the shell analysis, we now consider the inde- 
pendent dimension of metallicity (using color as a proxy) . 
Figure [14] shows that the stream GCs may have a differ- 
ent color distribution from the rest of the GCs. The 
former are narrowly confined to a blue color range of 
[g - i)o = 0.71-0.83 (a range in [Fe/H] of - 0.5 dex), 
while the latter include a similar blue peak but also a 
second, much redder peak, with [g — i)o = 0.97-1.13. 
From the K-S test, this difference has a 22% chance of 
occurring randomly. As with the shell, the color infor- 
mation of the stream provides independent confirmation 
of its distinctiveness. 

As discussed for the shell region, we could in principle 
connect the peak color of the stream GCs to a possible 
progenitor, but the statistics are currently too poor for 
strong constraints. Even so, we can still mention one ex- 
ample of a potential analogue to the stream progenitor: 
the Virgo dwarf elliptical IC 3328, with a stellar disper- 
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Fig. 13. — Globular clusters near outer "stream A", (a): positional spaee (a zoom-in of Figure[TIa)). The five candidate stream GCs are 
outlined in purple (the two easternmost ones are almost overlapping). The bar with arrows again shows a 100 kpe scale, (b): phase-space 
(using distance from the center of M87). Velocity uncertainties are indicated by error bars. One GC at upper right is probably bound to 
the low-luminosity SO NGC 4461 (red cross). The systemic velocity is shown by a horizontal line, and a diagonal line illustrates our best- fit 
stream model. 
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Fig. 14. — Color distribution of GCs in the stream A region. The 
two colored histograms show the stream and non-stream objects, 
per the leg end , with alternative smoothed distributions also shown 
(see Figure lTOf b)). The stream GCs do not include any red objects, 
but the difference with the non-stream GCs is not statistically sig- 
nificant. 



sion of 35 km s~^, hosting ^ 40 GCs with a color peak 
near {g — z)o ^ 0.92, i.e., {g — i)o ^ 0.77. 

Could this apparent GC stream from a dwarf galaxy be 
associated with the previously identified stellar stream? 
The latter has a luminosity of log [Ly / Lv,q) ~ 8.4 that 
is consistent with this hypothesis (compare the Ly-cTc 
trends for ACSVCS galaxies in the lower-left panel of 
Figure fTTj) . Its col or of (B — V) ~ .8 also suggests 
a dwarf progenitor (jRudick et al.l 120101 ) . although more 
massive progenitors are also possible within the uncer- 
tainties. 

A puzzle is presented by the spatial offset of 10 kpc 
between the stellar and GC substructures (Figure [TST b)). 
One possibility that will require further modeling is that 
the GCs and stars followed different trajectories owing 
to their different initial binding energies. The nucleus of 
the dwarf could also very well still be visible around M87, 
although we have not identified an obvious candidate. 

Another interesting possibility is that both the stream 
and the shell are part of the same accretion event, since 
their peak GC colors and velocity dispersions are simi- 
lar. Exploring this scenario will require more extensive 
simulations than we can attempt here. 

5. THEORETICAL ANALYSIS 

We now attempt to model the dynamics of the ob- 
served substructures, based mostly on numerical simula- 
tions but also making use of analytic representations of 
orbit kinematics. 

The narrow diagonal tracks of the stream and shell in 
phase space can be understood in simple terms as the 
near-radial infall of objects with similar initial potential 
energies. Analog ous features are fou nd in Local Group 
accretion events (jCilber t et al."2007*) and in simulations 
of minor mergers (e.g., iQuinn.1984; .Dupraz fc Combei 
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as follows. More massive satellites produce more diffuse 
streams owing to the higher internal velocity dispersions 
and larger sizes of these systems. In cuspy potentials, 
low angular momentum orbits can bring the satellites 
close to the potential center and lead to broader shells 
or fan-like features, while in orbits farther out, the 
weaker tidal shear produces narrower loops and streams. 
Finally, any dynamically cold substructure will tend to 
diffuse over time due either to nonaxisymmetry in the 
host potential or to perturbations from other satellites. 

As a proof of principle, we carry out a basic set of simu- 
lations inspired by M87 and its substructures. These are 
not intended to explore the (very large) orbital and pro- 
genitor parameter space for the infalling satellite galaxy, 
or to provide a unique, robust match to the data, but 
to illustrate qualitatively how comparable features can 
emerge in phase-space from the tidal stripping of infalling 
galaxies. 

We conduct TV-body simulations of self-gravitating low 
mass galaxies orbiting in a fixed spherical gravitational 
potential. For this cluster-mass potential, we use a 
Navarro-Frenk- White model for the da rk matter halo 
jNavarro et al.' '1996'; 'McLa ughlinl Il999f ). and embed a 
iHernauist (1990) profile to represent the stellar mass of 
the central galaxy. The relevant parameters are reported 
in Tabled 

For the inner shell, the main features that we wish to 
reproduce are the apparent chevron morphology in phase 
space including a low velocity dispersion, with a high 
degree of azimuthal mixing in real-space (see Figure [1]). 
For the outer stream, the goal is to find a thin structure in 
both real and phase-space, which is not closely connected 
to the progenitor galaxy or nucleus. 

For both features, we adopt the same low mass galaxy 
model for the infalling satellite, and consider a modest 
number of orbital configurations. The infalling satellite 
galaxy follows a single-component Hernquist model, with 
characteristic parameters summarized in Table [3] A real 
galaxy will consist of multiple components (stars, GCs, 
gas, and dark matter), but this introduces additional de- 
grees of freedom to the problem which we arc for now 
keeping as conceptually simple as possible. 

We place the satellite on elongated orbits around the 
potential center F^ and run the s imulations using an A^- 
body treecode (jHernqui"stl 119871 ). which simulates the 
fully self-gravitating response of the infalling satellite. 
The simulations were evolved for several Gyr using a 
fixed timestep of 0.25 Myr. We create real-space and 
phase-space plots of all the satellite particles for each 
time-step, using a small selection of viewing angles. We 
then select by eye the best combination of timestep and 
viewing angles that qualitatively resembles the data. 

Next, we sample a subset of the particles to represent 

Our adopted apocentric distances of '~ 100-200 kpc may not 
seem consistent with satellites that fall in from the large-scale envi- 
ronment, but our goal is to set up an idealized model of accretion 
dynamics once an obje ct ha s migrated inwards by sca ttering or 
dynamical friction (e.s.. iFaTtenbac her fc Mathew3l2005f l. 



the GC observations, plotting 10'* particles to obtain a 
general sense of the substructure morphology, and an ad- 
ditional subsample of 100 particles for more direct com- 
parison with the data. Our default selection consists of 
an unbiased, random sample, which can be considered 
as representing the limiting case where most of the dark 
halo has been stripped away, with the remaining self- 
gravitating material traced well by the stars and GCs. 
In almost every real galaxy, the metal-poor GCs occupy 
a more extended distribution than the field stars, so we 
have tried the alternative limiting case where the least 
bound particles are selected as "GCs" . The results from 
the two sampling schemes turn out to be qualitatively 
similar. 

Clearly, our procedures are only the initial steps to- 
ward modeling the M87 halo kinematics. Ideally, one 
would like to: use a live, triaxial potential that in- 
cludes the effects of dynamical friction and other sub- 
structures; explore a wide range of satellite and orbital 
parameters; and use an objective, quantitative metric 
for optimizing the model fits. However, our current ap- 
proach is adequate for providing a broad sense of some of 
the types of substructure signatures that are physically 
plausible , which one can see in related analyses of other 
syste ms (IHern auist fc QuiimI 19881: iMerrifield fc KuiikenI 
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Gilbert et al.ii2007l: iFardal et al.ll2007t iPenarrubia et al l 
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The best qualitative matches that we identified for the 
shell and the stream (separately) were shown in Figure[2j 
The shell progenitor has made a close passage to the 
center of M87 and was largely disrupted, leaving a clas- 
sic chevron signature in phase space, while also being 
spread out in positional space. The stream progenitor 
has a larger pericenter and is disrupted more gradually, 
maintaining integrity as a narrow stream in both phase 
space and positional space for many orbits. 

In a real system, any triaxiality of the potential would 
cause the orbits to phase- mix more rapidly, and the time- 
evolution of the mass distribution would presumably 
smear out the cold substructures even further (as well as 
inhibit any orbital resonances that might in principle be 
an alternative mechanism for producing cold kinematical 
features). Simulations of cluster-sca le tidal streams in a 
cosmological context (jRudick et al.l 12009^ indicate that 
streams are disrupted over timescales of ~ 2-4 times the 
local dynamical timescale, idyn ~ t/'^c- For the shell and 
stream, t^yn 0.2 and 0.5 Gyr, respectively, so their 
most likely lifetimes are '-^ 0.5 Gyr and ~ 1.5 Gyr. 

One important feature of coherent substructures is 
that they trace out a close family of orbits over large 
scales in galaxy halos, and can thus be used to pro- 
vide unique constraints on the gravitational potential, 
including the radi al distribution and shape of the dark 
matter halo fe.g..lLaw et al.l l2009l: iKoposov et all 120101: 
iVarghese et allboilD . The case of a spherical system is 
conceptually simple: the potential energy lost by parti- 
cles as they fall in towards the center of the host galaxy 
is reflected in a monotonic trend of increasing veloc- 
ities with decreasing radiu s, modulo projection effects 
(|Merrifleld fc Kuiikenl[T99l . 

We derive a new formulation of this method, assuming 
a spherical gravitational field of a power-law form, where 
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TABLE 3 

Parameters of Accretion Simulations 



Cold shell Cold stream Hot shell 



satellite effective radius, i?e [kpc] 


1 


1 


4 


satellite velocity dispersion, ap(Re) [km s~^] 


35 


35 


200 


satellite mass, Msat [A^q] 


1.8 X 10^ 


1.8 X 10^ 


3.2 X 10" 


particle number, Np 


65,536 


65,536 


65,536 


particle mass, Mp [Mq] 


2.75 X 10"* 


2.75 X 10* 


4.9 X 10*^ 


gravitational softening length, rj [pc] 


20 


20 


22 


initial transverse velocity, vt / Vc 


0.05 


0.1 


0.1 


initial apocentric radius, rapo [kpc] 


90 


200 


90 


initial pericentric radius, rperi [kpc] 


2 


12 


4.7 


host dark halo virial mass, M200 [Af©] 


4.2 X lO^"' 


4.2 X 10^"' 


4.2 X 10" 


host dark halo virial radius, r2oo [Mpc] 


1.55 


1.55 


1.55 


host dark halo scale radius, [Mpc] 


0.564 


0.564 


0.564 


host stellar effective radius, -Re [kpc] 


7 


7 


7 


host circular velocity at i?c, Vc [km s~^] 


475 


475 


475 


timestep, [Myr] 


0.23 


0.23 


0.14 



the total density, circular velocity, and potential are: 
Ptot oc r"", 



Jo 



$(r) 



(9) 
(10) 

(11) 



where a 
100 kpc 



< 2, as is probably appropriate on ~ 10- 

scales in M87 (cf. also section 2.2(f) of 

Binnev fc T remainc 2008). Following equation (3) of 



Merrifield fc Kuiiken (1998i) . we consider a thin shell of 



radius rapo where particles have zero velocity (i.e., they 
are on radial orbits at apocenter). By conservation of 
energy and by geometry, their line-of-sight velocities at 
3-D radius r and 2-D radius R are then given by 



vlos =21 



n2 

^ ) 



ap 



(12) 



The velocity "edge" of the shell is then determined by 
solving for the maximum of fLos a-t a given radius R. 
By fitting this simple model to a cold phase-space fea- 
ture, we can then estimate the value of Vc near rapo- An 
even more simplified way to understand this constraint 
(jMerrifield fc Kuijkcn 1998) is that the slope of the shell 
velocity edge with radius in phase-space is approximately 
equal to fc/rapo- In practice, the dominant effect is the 
change of potential energy, so the determination of Vc is 
fairly insensitive to the form adopted for $(r). 

In Figure [15] we show this model applied to the 
stream and to two possible shells in M87, where we have 
taken the maximum-likelihood linear-slope fits from Sec- 
tions 13.41 and H] as starting points for approximate fits 
in phase-space. The model reproduces nicely the ob- 
served radial trends, and suggests some possible con- 
nections in phase-space between the shells and nearby 
low-luminosity ellipticals. Also shown in the figure are 
halo PNe, which show hints of tracing some of the 
same substructures seen in the GCs. IVIore precise or- 
bital fits to the substructures and the potential progen- 
itors will require d edicated dynamical modeling _(e.g., 
iRomanowskv fc Kochanck.2001; .Murphv et al...201l!) . 
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Fig. 15. — Phase-space diagram of M87 as in Figure [Jib), with 
three monoenergetic model curves over-plotted (see text for de- 
tails). GC and satellite galaxy data are show n by orange dots and 
red crosses as before, with planetary nebulae IjDohertv et al."2009l) 
also now included as green squares. The outer stream subgroup of 
GCs (at Rm ^ 35 arcmin) is highlighted with larger circles. The 
inner shell and outer stream GCs can be represented well by model 
curves as shown, and there is a hint of an intermediate-radius shell 
as marked, whose apex could explain the low velocity dispersion 
found in a subgroup of the PNe. 
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Fig. 16. — Cold accretion model with higher initial angular momentum, viewed near face-on (compare Figure [J] note that both the 
vertical and horizontal scales have been shrunk in tandem for a fair comparison of angles). The snapshot corresponds to 2 Gyr after the 
initial apocenter. 

The opening angle of the chevron is smaller than in Figure [^b), and in position-space the stream stays highly coherent for many orbits 

(in the absence of perturbations in the potential). 



There is a hitch with these model curves: the imphed 
circular velocities differ radically from typical estimates 
for the M87 potential. At r - 20, 32, 41 arcmin 90, 
150, 200 kpc), we infer ~ 270, 1200, 1400 km s'^ 
albeit with substantial uncertainties. Previous estimates 
from dynamical and X-ray analyses yielded fc 650- 
900 km s~^ over this e ntire radial range ()Das et al.l 
[20T0I: IMurohv et aPIMTh . while our own equilibrium- 
based dynamical analysis indicates Vc ~ 400-550 km s~^ 
(S-|-ll)l3 These discrepancies should not be a major 
concern at this point because there are various complica- 
tions to consider including the radius conv ention of Sec- 
tion 13.11 as well as non-spherical effects (jjflkova et al.l 
[2OIQ1) . In particular, a planar orbit viewed near face-on 
will show depressed line-of-sight velocities, and a reduced 
opening angle of the chevron. 

An illustration of the latter effect is provided in Fig- 
ure [161 which shows an alternative orbit for the shell that 
has higher angular momentum than our fiducial case, and 
maintains better planar integrity. Here the "measured" 
Vc would be ~ 250 km s~^, despite the true value in the 
model being 635 kms~^. Note that the simulation shown 
is also representative of a large family of orbits that are 
not a good match to the shell observations because of the 
high degree of spatial coherence. 

As discussed at length in Section [5?^ there is an appar- 
ent conflict between the shell's inferred kinematical cold- 
ness and the large number of GCs, since these imply very 
different progenitor masses. We consider this issue more 

The local escape velocity is at least We regard 

the Vc ~ 550 km model as plausibly describing a massive 
group-sized dark matter halo around M87, in which case Ve > 
800 km . This would indicate a bound- velocity range of ~ 500— 
2100 km s^^, and probably wider, suggesting that most or all of 
the objects in Figure [Tsl are bound to M87 (modulo the unknown 
proper-motion velocities). 



explicitly by changing the simulated shell progenitor to 
correspond to a more massive galaxy, with parameters 
summarized in Table |31 

We show an example post-disruption snapshot in Fig- 
ure 1171 finding as expected that the resulting shell kine- 
matics are hotter than in the dwarf-galaxy simulation, 
and do not compare as favorably to the observations 
(compare Figures [2{b) and [31[b)). However, as previ- 
ously discussed, this tension is alleviated if we modify 
our default radius convention and chevron-based model 
(compare Figure [3l[a)). 

These simulations bracket the simple cases of very low- 
and high-mass progenitor galaxies, but a worthwhile en- 
deavor for the future would be simulating more detailed 
intermediate-mass cases, including bulge, disk, stellar 
halo, and dark matter halo subcomponents. An addi- 
tional aspect that could be pursued is the modeling of 
velocity gradients with azimuthal angle, which should 
be ge nerically found in tidal features (e.g., Helm i et al.l 
120031 ). Our initial look at the M87 shell data shows no 
clear sign of such effects. 

6. SUMMARY 

We have used an unprecedentedly wide-field data set 
of high-precision globular cluster velocities to search 
for phase-space substructure around the central cluster 
galaxy M87. The large majority of the data are from 
S-l-11, with a few additional velocities obtained from the 
literature, and from a new sample of Keck/DEIMOS 
spectroscopy presented here that targets a known fila- 
ment of stellar light in the outer halo. Analysis oi (g — i) 
colors and spectroscopic metallicities for a subset of our 
GC sample confirms that the colors provide a good proxy 
for metallicity, allowing for chemical tagging of substruc- 
tures. 

We identify two candidate features: one a small stream 
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Fig. 17. — Hot accretion model, similar to the shell simulation in Figure[2]but with a lOOx more massive progenitor. Here the snapshot 
corresponds to 1.5 Gyr after the initial apocenter. Similar features are present as with the low-mass galaxy but are less distinct. 



near to the outer- halo stellar filament, and the other a 
previously unknown colossal "shell" in the inner halo. 
We have carried out extensive statistical tests and char- 
acterizations of these substructures. 

The outer GC stream appears to be a collection of 
15 ± 10 objects (after extrapolating to fainter magni- 
tudes) residing in a ~ 5 x 20 kpc clump with an in- 
ternal velocity dispersion of ^ 30 ± 20 km s^^. We 
have found the presence of this stream to be significant 
both in a group-finding analysis (p = 0.09), from a max- 
imum likelihood model {p — 0.025), and from its distinc- 
tive color distribution compared to the surrounding GCs 
(p = 0.22). 

The stream dispersion, its GC colors and numbers 
counts, the color of the stellar filament, and the fila- 
ment's narrowness all suggest the accretion of a dwarf 
galaxy. Low-mass kinematical substructures in elliptical 
galaxies have been inferred in other studies (jCote et al.l 
l200a IWoodlev fc Ha^[20Tltl . but in this case the as- 
sociation with a visible counterpart should allow for a 
much clearer determination of the accretion dynamics. 

The shell of GCs has a chevron-like shape in phase- 
space that rese mbles classic expectations f or a disrupted 
infalling system (|Hernquist fc Qu inn"'1988^ , although the 
interpretation in detail is more challenging. The appear- 
ance of the shell is sharpened by a change of coordi- 
nates from simple radius to elliptical, circular-equivalent 
radius. The significance of the shell in phase-space is 
established using a group-finding algorithm {p < 0.01), 
a crude entropy metric {p = 0.31), and a maximum- 
likelihood chevron-based model {p < 0.01 compared to 
a pure Gaussian model). The color distribution of the 
shell and non-shell GCs from the same region is differ- 
ent {p = 0.13), providing independent evidence for the 
existence of the shell. The phase-space morphology re- 
sembles a simple chevron model more than a wedge or 
tapered Gaussian, but may very well have a complex, in- 
termediate morphology; however, such uncertainties do 



not compromise our basic finding of a relatively massive 
and cold substructure. 

We find from maximum-likelihood fitting that the shell 
contains 500-1100 GCs, which is ~ 5 times the entire 
Milky Way system, and comprises around 20% of M87's 
GCs in the 50-95 kpc region. This number suggests a 
stellar debris field of {2-6)x10^'^Lq, brought in by ei- 
ther a large group of dwarf galaxies, or by a single giant 
elliptical or lenticular galaxy. 

The dramatic appearance of the shell in phase-space 
is due to its relatively cold kinematics, with an esti- 
mated velocity dispersion of 10-30 km (or up to 
^ 120 km if we make generous allowances for uncer- 
tainties concerning the shell morphology). The lower ve- 
locity dispersion estimates are difficult to reconcile with 
the number of inferred GCs, as they most naturally imply 
a dwarf galaxy hosting only a handful of GCs. However, 
after allowing for systematic uncertainties in our analy- 
ses, we suggest a solution involving an E/SO progenitor 
with a luminosity of ^ 0.5 L*, which could very well also 
be the parent galaxy for one of the tidally-stripped galax- 
ies now visible close in to M87. 

There is a critical need for additional data from the 
shell region, including new velocities of both GCs and 
PNe, in order to confirm the presence of the substruc- 
ture and to refine our estimates of its characteristics. 
Along these lines, we have recently obtained another cy- 
cle of GC spectroscopy from MMT/Hectospec that in- 
cludes ~ 35 velocity measurements of independent ob- 
jects in the shell region. A preliminary analysis yields no 
obvious sign of substructure, but the statistics from the 
new data are poor, and at least another ^ 100 velocities 
are needed. 

We have carried out simplified iV-body simulations of 
satellite galaxy accretion in a static cluster -l-BCG po- 
tential, finding that these readily produce streams and 
chevrons in phase space. In one simulation, a dwarf 
galaxy has recently fallen in and is partially disrupted. 
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forming a classic, narrow tidal stream in both positional 
and phase space. In a second case, the initial infall orbit 
had much lower angular momentum, so that the satellite 
galaxy passed close to the center of M87 and was quickly 
disrupted. Its remnants are strewn more haphazardly 
across positional space and could be difBcult to discern 
against the main body of M87, while a distinct chevron 
persists in phase space, with the apex corresponding to 
a shell turning point. 

These simulations compare favorably to the M87 ob- 
servations because the low internal velocity dispersions 
of the dwarf progenitors result in cold kinematics for the 
tidal debris. Simulations of a more massive progenitor 
produce hotter kinematics that do not as strongly re- 
semble the data for the "shell", reiterating the tension 
between the different constraints on the shell progenitor 
mass, which remains unresolved pending further observa- 
tions and modeling. It would be worth also investigating 
alternative explanations for the substructure, such as a 
shell-like boundary between well-mixed GC subpopula- 
tions formed or accreted in different epochs. 

We have explored the use of the velocity-radius slopes 
of the substructures in phase-space estimate the underly- 
ing gravitational potential and dark matter distribution 
around MS 7. However, the interpretation of the data 
appears to be complicated by non-spherical effects. 

Taken at face value, the shell GCs imply a substan- 
tial, recent accretion event in the halo of M87. This is 
the first large substructure identified around a central 
cluster galaxy, the first large stellar substructure with a 
clear kinematical detection in any type of galaxy beyond 
the Local Group, and by far the largest shell or stream 
discovered in the Local Universe. 

M87 was a favored target for our survey as well as 
previous ones because of its proximity, rich GC system, 
and general lack of obvious dynamical disturbance. The 
large shell we have discovered is an example of the lurking 
accretion signatures that could be found from detailed 
phase-space studies of apparently placid galaxies. 

There are however more indirect lines of evidence for a 
recent minor merger in M 87, such as an off-ce nter su- 
permassive black hole f)Batcheldor et al.l 120101) , a lop- 
sided central GC velocity distribution (S-l-11), an accre- 
tion powered jet, and several surrounding peculiar low- 
luminosity satellites that could be the nuclei of since- 
shredded disk galaxies. Further theoretical work is 
needed to see if these features could all be explained by a 
unified scenario involving a recent gas-rich minor merger. 

The shell and stream in M87 have probable lifetimes 
of ^ 0.5 Gyr and ^ 1.5 Gyr, respectively. Given our 
best guesses for the progenitors of these features, the 
total luminosity of M87 would be built up by ~ 10 or 
~ 1000 events corresponding to the shell or the stream, 
respectively. If these recent events are representative of 
the long-term evolution of M87, and by extension of the 
class of brightest cluster galaxies in general, then they 
support a picture of late-epoch gro wth by hierarchical 
assembly ()De Lucia &: Blaizotir2007f) . 

More generally, the M87 results are another demon- 
stration of the potency of GCs for providing unique in- 
formation about the stellar halos of galaxies. We note 



Arnold, J. A., Romanowsky, A. J., Brodie, J. P., Chomiuk, L 
Spitler, L. R., Strader, J., Benson, A. J., & Forbes, D. A. 2011, 
ApJ, 736, L26 



the concluding prediction of iNaab et al.l ()2009D for their 
landmark two-phase galaxy formation scenario: "...the 
outer parts of massive giant ellipticals will tend to be 
old, blue, metal-poor, and relatively uniform from galaxy 
to galaxy since they are all composed essentially of the 
debris from tidally destroyed accreted small systems." 
Globular clusters have in fact long provided evidence for 
exactly this scenario. 
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